Changeset 21381


Ignore:
Timestamp:
11/16/16 18:58:41 (8 years ago)
Author:
felicity
Message:

CHG: ESTAR in terms of epseff and B

Location:
issm/trunk-jpl/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r21265 r21381  
    4848        IssmDouble vx,vy,vz,vmag;
    4949        IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3];
    50         IssmDouble epso,epsprime;
     50        IssmDouble epseff,epsprime;
    5151        int         dim;
    5252        IssmDouble *xyz_list = NULL;
    … …  
    101101                        dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
    102102                }
    103                 EstarStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
    104                 lambdas[iv]=EstarLambdaS(epso,epsprime);
     103                EstarStrainrateQuantities(&epseff,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
     104                lambdas[iv]=EstarLambdaS(epseff,epsprime);
    105105        }
    106106
    … …  
    15141514                                name==MaterialsRheologyBbarEnum ||
    15151515                                name==MaterialsRheologyNEnum ||
    1516                                 name==MaterialsRheologyKoEnum ||
    1517                                 name==MaterialsRheologyKobarEnum ||
    15181516                                name==MaterialsRheologyEcEnum ||
    15191517                                name==MaterialsRheologyEcbarEnum ||

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r21206 r21381  
    22892289                        break;
    22902290                case MatestarEnum:
    2291                         this->InputDepthAverageAtBase(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum);
     2291                        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    22922292                        this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
    22932293                        this->InputDepthAverageAtBase(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum);

  • issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp

    r21112 r21381  
    134134/*}}}*/
    135135IssmDouble Matestar::GetA(){/*{{{*/
    136         _error_("not implemented yet");
     136        /*
     137         * A = 1/B^n
     138    */
     139
     140        IssmDouble B,n;
     141
     142        element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
     143        n=this->GetN();
     144
     145        return pow(B,-n);
    137146}
    138147/*}}}*/
    139148IssmDouble Matestar::GetAbar(){/*{{{*/
    140         _error_("not implemented yet");
     149        /*
     150         * A = 1/B^n
     151         */
     152
     153        IssmDouble B,n;
     154
     155        element->inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
     156        n=this->GetN();
     157
     158        return pow(B,-n);
    141159}
    142160/*}}}*/
    143161IssmDouble Matestar::GetB(){/*{{{*/
    144         _error_("not implemented yet");
     162        /*Output*/
     163        IssmDouble B;
     164
     165        element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
     166        return B;
    145167}
    146168/*}}}*/
    147169IssmDouble Matestar::GetBbar(){/*{{{*/
    148 
    149         _error_("not implemented yet");
     170        /*Output*/
     171        IssmDouble Bbar;
     172
     173        element->inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum);
     174        return Bbar;
    150175}
    151176/*}}}*/
    … …  
    159184}
    160185/*}}}*/
     186IssmDouble Matestar::GetEc(){/*{{{*/
     187
     188        /*Output*/
     189        IssmDouble Ec;
     190
     191        element->inputs->GetInputAverage(&Ec,MaterialsRheologyEcEnum);
     192        return Ec;
     193}
     194/*}}}*/
     195IssmDouble Matestar::GetEs(){/*{{{*/
     196
     197        /*Output*/
     198        IssmDouble Es;
     199
     200        element->inputs->GetInputAverage(&Es,MaterialsRheologyEsEnum);
     201        return Es;
     202}
     203/*}}}*/
    161204IssmDouble Matestar::GetN(){/*{{{*/
    162         _error_("not implemented yet");
     205        return 3.;
    163206}
    164207/*}}}*/
    … …  
    228271}
    229272/*}}}*/
    230 IssmDouble Matestar::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
     273IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
    231274
    232275        /*Intermediaries*/
    233276        IssmDouble viscosity;
    234         IssmDouble E,lambdas;
    235         IssmDouble epso,epsprime_norm;
     277        IssmDouble epseff,epsprime_norm;
     278        IssmDouble lambdas;
    236279        IssmDouble vmag,dvmag[3];
     280        IssmDouble B,Ec,Es,E;
    237281
    238282        /*Calculate velocity magnitude and its derivative*/
    … …  
    249293        }
    250294
    251         EstarStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
    252         lambdas=EstarLambdaS(epso,epsprime_norm);
     295        EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
     296        lambdas=EstarLambdaS(epseff,epsprime_norm);
     297
     298        /*Get B and enhancement*/
     299        B=GetB(); _assert_(B>0.);
     300        Ec=GetEc(); _assert_(Ec>=0.);
     301        Es=GetEs(); _assert_(Es>=0.);
    253302
    254303        /*Get total enhancement factor E(lambdas)*/
    255         E = Ec + (Es-Ec)*lambdas*lambdas;
     304        E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
    256305
    257306        /*Compute viscosity*/
    258         _assert_(E>0.);
    259         _assert_(ko>0.);
    260         _assert_(epso>0.);
    261         viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));
     307        /*if no strain rate, return maximum viscosity*/
     308        if(epseff==0.){
     309                viscosity = 1.e+14/2.;
     310                //viscosity = B;
     311                //viscosity=2.5*pow(10.,17);
     312        }
     313        else{
     314                viscosity = B/(2.*pow(E*epseff*epseff,1./3.));
     315        }
    262316
    263317        /*Assign output pointer*/
    … …  
    324378        IssmDouble vx,vy,vz;
    325379        IssmDouble dvx[3],dvy[3],dvz[3];
    326         IssmDouble ko,Ec,Es;
     380        IssmDouble B,Ec,Es;
    327381
    328382        /*Get velocity derivatives in all directions*/
    … …  
    344398        }
    345399
    346         /*Get enhancement factors and ko*/
    347         Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
    348         Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
    349         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
    350         ec_input->GetInputValue(&Ec,gauss);
    351         es_input->GetInputValue(&Es,gauss);
    352         ko_input->GetInputValue(&ko,gauss);
    353 
    354400        /*Compute viscosity*/
    355         *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     401        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
    356402}
    357403/*}}}*/
    … …  
    364410        IssmDouble vx,vy,vz;
    365411        IssmDouble dvx[3],dvy[3],dvz[3];
    366         IssmDouble ko,Ec,Es;
    367412
    368413        /*Get velocity derivatives in all directions*/
    … …  
    384429        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
    385430
    386         /*Get enhancement factors and ko*/
    387         Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
    388         Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
    389         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
    390         ec_input->GetInputValue(&Ec,gauss);
    391         es_input->GetInputValue(&Es,gauss);
    392         ko_input->GetInputValue(&ko,gauss);
    393 
    394431        /*Compute viscosity*/
    395         *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     432        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
    396433}/*}}}*/
    397434void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    … …  
    405442        IssmDouble vx,vy,vz;
    406443        IssmDouble dvx[3],dvy[3],dvz[3];
    407         IssmDouble ko,Ec,Es;
     444        IssmDouble B,Ec,Es;
    408445
    409446        /*Get velocity derivatives in all directions*/
    … …  
    428465        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
    429466
    430         /*Get enhancement factors and ko*/
    431         Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
    432         Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
    433         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
    434         ec_input->GetInputValue(&Ec,gauss);
    435         es_input->GetInputValue(&Es,gauss);
    436         ko_input->GetInputValue(&ko,gauss);
    437 
    438467        /*Compute viscosity*/
    439         *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     468        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
    440469}/*}}}*/
    441470void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/

  • issm/trunk-jpl/src/c/classes/Materials/Matestar.h

    r20945 r21381  
    6969                IssmDouble GetD();
    7070                IssmDouble GetDbar();
     71                IssmDouble GetEc();
     72                IssmDouble GetEs();
    7173                IssmDouble GetN();
    7274                bool       IsDamage();
    … …  
    8284                void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    8385                /*}}}*/
    84                 IssmDouble GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
     86                IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
    8587};
    8688

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r20690 r21381  
    8080                        break;
    8181                case MatestarEnum:
    82                         iomodel->FetchDataToInput(elements,"md.materials.rheology_ko",MaterialsRheologyKoEnum);
     82                        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
    8383                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
    8484                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
    … …  
    8686                        switch(iomodel->domaindim){
    8787                                case 2:
    88                                         elements->InputDuplicate(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum);
     88                                        elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    8989                                        elements->InputDuplicate(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
    9090                                        elements->InputDuplicate(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum);

  • issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp

    r20945 r21381  
    33#include "../Exceptions/exceptions.h"
    44#include "./elements.h"
    5 void EstarStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
     5void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
    66
    77        /*Intermediaries*/
    88        IssmDouble omega[3];
    99        IssmDouble nrsp[3],nrsp_norm;
    10         IssmDouble eps[3][3],epso;
     10        IssmDouble eps[3][3],epseff;
    1111        IssmDouble epsprime[3],epsprime_norm;
    1212
    … …  
    2424                nrsp[0] = 0.;
    2525                nrsp[1] = 0.;
    26                 nrsp[2] = 1.;
     26                nrsp[2] = 0.;
    2727        }
    2828        else{
    … …  
    3737        eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];
    3838
    39         /*octahedral strain rate*/
    40         epso = 0.;
    41         for(int i=0;i<3;i++) for(int j=0;j<3;j++) epso += eps[i][j]*eps[i][j];
    42         if(epso==0.) epso=1.e-14;
    43         epso=sqrt(epso)/sqrt(3.);
     39        /*effective strain rate*/
     40        epseff = 0.;
     41        /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
     42        epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] +  eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]);
    4443
    45         /*Compute the shear strain rate on the non ratating shear plane*/
     44        /*Compute the shear strain rate on the non rotating shear plane*/
    4645        epsprime[0]=0.;
    4746        epsprime[1]=0.;
    … …  
    7473       
    7574        /*Assign output pointers*/
    76         *pepso=epso;
     75        *pepseff=epseff;
    7776        *pepsprime_norm=epsprime_norm;
    7877}/*}}}*/
    … …  
    114113                omega[0] = 0.;
    115114                omega[1] = 0.;
    116                 omega[2] = 1.;
     115                omega[2] = 0.;
    117116        }
    118117        else{
    … …  
    123122
    124123}/*}}}*/
    125 IssmDouble EstarLambdaS(IssmDouble epso, IssmDouble epsprime_norm){/*{{{*/
    126    _assert_(epso>0.);
     124IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm){/*{{{*/
     125        IssmDouble lambdas;
     126
    127127        _assert_(epsprime_norm>=0.);
    128         return sqrt(2*epsprime_norm*epsprime_norm/(3*epso*epso));
     128        if(epseff==0.){
     129                lambdas=0.;
     130        }
     131        else{
     132                lambdas=sqrt(epsprime_norm*epsprime_norm/(epseff*epseff));
     133        }
     134        return lambdas;
    129135}/*}}}*/

  • issm/trunk-jpl/src/c/shared/Elements/elements.h

    r20687 r21381  
    1515IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat);
    1616// IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
    17 IssmDouble EstarLambdaS(IssmDouble epso, IssmDouble epsprime_norm);
     17IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm);
    1818void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag);
    19 void EstarStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
     19void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
    2020IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
    2121                                 IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm,

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r21344 r21381  
    214214        MaterialsRheologyLawEnum,
    215215        MaterialsRheologyNEnum,
    216         MaterialsRheologyKoEnum,
    217         MaterialsRheologyKobarEnum,
    218216        MaterialsRheologyEcEnum,
    219217        MaterialsRheologyEcbarEnum,

  • issm/trunk-jpl/src/m/classes/matestar.m

    r21377 r21381  
    1818                mixed_layer_capacity       = 0.;
    1919                thermal_exchange_velocity  = 0.;
    20                 rheology_ko   = NaN;
     20                rheology_B    = NaN;
    2121                rheology_Ec   = NaN;
    2222                rheology_Es   = NaN;
    … …  
    3535        methods
    3636                function self = extrude(self,md) % {{{
    37                         self.rheology_ko=project3d(md,'vector',self.rheology_ko,'type','node');
     37                        self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
    3838                        self.rheology_Ec=project3d(md,'vector',self.rheology_Ec,'type','node');
    3939                        self.rheology_Es=project3d(md,'vector',self.rheology_Es,'type','node');
    … …  
    5858                                                self.rheology_Es = ones(size(inputstruct.rheology_B));
    5959                                                self.rheology_Ec = ones(size(inputstruct.rheology_B));
    60                                                 self.rheology_ko = inputstruct.rheology_B.^(-3)*3/2;
     60                                                self.rheology_B = inputstruct.rheology_B;
    6161                                        end
    6262                                otherwise
    … …  
    121121                        md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    122122                        md = checkfield(md,'fieldname','materials.mu_water','>',0);
    123                         md = checkfield(md,'fieldname','materials.rheology_ko','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
     123                        md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    124124                        md = checkfield(md,'fieldname','materials.rheology_Ec','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    125125                        md = checkfield(md,'fieldname','materials.rheology_Es','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    … …  
    152152                        fielddisplay(self,'mixed_layer_capacity','mixed layer capacity [W/kg/K]');
    153153                        fielddisplay(self,'thermal_exchange_velocity','thermal exchange velocity [m/s]');
    154                         fielddisplay(self,'rheology_ko','octahedral flow law parameter [s^-1 Pa^-2]');
     154                        fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/3)]');
    155155                        fielddisplay(self,'rheology_Ec','compressive enhancement factor');
    156156                        fielddisplay(self,'rheology_Es','shear enhancement factor');
    … …  
    176176                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double');
    177177                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
    178                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_ko','format','DoubleMat','mattype',1);
     178                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1);
    179179                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Ec','format','DoubleMat','mattype',1);
    180180                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Es','format','DoubleMat','mattype',1);
    … …  
    202202                        writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity);
    203203                        writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
    204                         writejs1Darray(fid,[modelname '.materials.rheology_ko'],self.rheology_ko);
     204                        writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B);
    205205                        writejs1Darray(fid,[modelname '.materials.rheology_Ec'],self.rheology_Ec);
    206206                        writejs1Darray(fid,[modelname '.materials.rheology_Es'],self.rheology_Es);
Note: See TracChangeset for help on using the changeset viewer.