Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21380) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21381) @@ -48,5 +48,5 @@ IssmDouble vx,vy,vz,vmag; IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3]; - IssmDouble epso,epsprime; + IssmDouble epseff,epsprime; int dim; IssmDouble *xyz_list = NULL; @@ -101,6 +101,6 @@ dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); } - EstarStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); - lambdas[iv]=EstarLambdaS(epso,epsprime); + EstarStrainrateQuantities(&epseff,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); + lambdas[iv]=EstarLambdaS(epseff,epsprime); } @@ -1514,6 +1514,4 @@ name==MaterialsRheologyBbarEnum || name==MaterialsRheologyNEnum || - name==MaterialsRheologyKoEnum || - name==MaterialsRheologyKobarEnum || name==MaterialsRheologyEcEnum || name==MaterialsRheologyEcbarEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21380) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21381) @@ -2289,5 +2289,5 @@ break; case MatestarEnum: - this->InputDepthAverageAtBase(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum); + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum); this->InputDepthAverageAtBase(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum); Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21380) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21381) @@ -134,18 +134,43 @@ /*}}}*/ IssmDouble Matestar::GetA(){/*{{{*/ - _error_("not implemented yet"); + /* + * A = 1/B^n + */ + + IssmDouble B,n; + + element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum); + n=this->GetN(); + + return pow(B,-n); } /*}}}*/ IssmDouble Matestar::GetAbar(){/*{{{*/ - _error_("not implemented yet"); + /* + * A = 1/B^n + */ + + IssmDouble B,n; + + element->inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum); + n=this->GetN(); + + return pow(B,-n); } /*}}}*/ IssmDouble Matestar::GetB(){/*{{{*/ - _error_("not implemented yet"); + /*Output*/ + IssmDouble B; + + element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum); + return B; } /*}}}*/ IssmDouble Matestar::GetBbar(){/*{{{*/ - - _error_("not implemented yet"); + /*Output*/ + IssmDouble Bbar; + + element->inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum); + return Bbar; } /*}}}*/ @@ -159,6 +184,24 @@ } /*}}}*/ +IssmDouble Matestar::GetEc(){/*{{{*/ + + /*Output*/ + IssmDouble Ec; + + element->inputs->GetInputAverage(&Ec,MaterialsRheologyEcEnum); + return Ec; +} +/*}}}*/ +IssmDouble Matestar::GetEs(){/*{{{*/ + + /*Output*/ + IssmDouble Es; + + element->inputs->GetInputAverage(&Es,MaterialsRheologyEsEnum); + return Es; +} +/*}}}*/ IssmDouble Matestar::GetN(){/*{{{*/ - _error_("not implemented yet"); + return 3.; } /*}}}*/ @@ -228,11 +271,12 @@ } /*}}}*/ -IssmDouble Matestar::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/ +IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/ /*Intermediaries*/ IssmDouble viscosity; - IssmDouble E,lambdas; - IssmDouble epso,epsprime_norm; + IssmDouble epseff,epsprime_norm; + IssmDouble lambdas; IssmDouble vmag,dvmag[3]; + IssmDouble B,Ec,Es,E; /*Calculate velocity magnitude and its derivative*/ @@ -249,15 +293,25 @@ } - EstarStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); - lambdas=EstarLambdaS(epso,epsprime_norm); + EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); + lambdas=EstarLambdaS(epseff,epsprime_norm); + + /*Get B and enhancement*/ + B=GetB(); _assert_(B>0.); + Ec=GetEc(); _assert_(Ec>=0.); + Es=GetEs(); _assert_(Es>=0.); /*Get total enhancement factor E(lambdas)*/ - E = Ec + (Es-Ec)*lambdas*lambdas; + E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.); /*Compute viscosity*/ - _assert_(E>0.); - _assert_(ko>0.); - _assert_(epso>0.); - viscosity = 1./(2*pow(ko*E*epso*epso,1./3.)); + /*if no strain rate, return maximum viscosity*/ + if(epseff==0.){ + viscosity = 1.e+14/2.; + //viscosity = B; + //viscosity=2.5*pow(10.,17); + } + else{ + viscosity = B/(2.*pow(E*epseff*epseff,1./3.)); + } /*Assign output pointer*/ @@ -324,5 +378,5 @@ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble ko,Ec,Es; + IssmDouble B,Ec,Es; /*Get velocity derivatives in all directions*/ @@ -344,14 +398,6 @@ } - /*Get enhancement factors and ko*/ - Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input); - Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input); - Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input); - ec_input->GetInputValue(&Ec,gauss); - es_input->GetInputValue(&Es,gauss); - ko_input->GetInputValue(&ko,gauss); - /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); } /*}}}*/ @@ -364,5 +410,4 @@ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble ko,Ec,Es; /*Get velocity derivatives in all directions*/ @@ -384,14 +429,6 @@ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; - /*Get enhancement factors and ko*/ - Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input); - Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input); - Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input); - ec_input->GetInputValue(&Ec,gauss); - es_input->GetInputValue(&Es,gauss); - ko_input->GetInputValue(&ko,gauss); - /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); }/*}}}*/ void Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ @@ -405,5 +442,5 @@ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble ko,Ec,Es; + IssmDouble B,Ec,Es; /*Get velocity derivatives in all directions*/ @@ -428,14 +465,6 @@ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; - /*Get enhancement factors and ko*/ - Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input); - Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input); - Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input); - ec_input->GetInputValue(&Ec,gauss); - es_input->GetInputValue(&Es,gauss); - ko_input->GetInputValue(&ko,gauss); - /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); }/*}}}*/ void Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21380) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21381) @@ -69,4 +69,6 @@ IssmDouble GetD(); IssmDouble GetDbar(); + IssmDouble GetEc(); + IssmDouble GetEs(); IssmDouble GetN(); bool IsDamage(); @@ -82,5 +84,5 @@ void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); /*}}}*/ - IssmDouble GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz); + IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz); }; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 21380) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 21381) @@ -80,5 +80,5 @@ break; case MatestarEnum: - iomodel->FetchDataToInput(elements,"md.materials.rheology_ko",MaterialsRheologyKoEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum); iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum); @@ -86,5 +86,5 @@ switch(iomodel->domaindim){ case 2: - elements->InputDuplicate(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum); + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); elements->InputDuplicate(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum); elements->InputDuplicate(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum); Index: /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21380) +++ /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21381) @@ -3,10 +3,10 @@ #include "../Exceptions/exceptions.h" #include "./elements.h" -void EstarStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/ +void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/ /*Intermediaries*/ IssmDouble omega[3]; IssmDouble nrsp[3],nrsp_norm; - IssmDouble eps[3][3],epso; + IssmDouble eps[3][3],epseff; IssmDouble epsprime[3],epsprime_norm; @@ -24,5 +24,5 @@ nrsp[0] = 0.; nrsp[1] = 0.; - nrsp[2] = 1.; + nrsp[2] = 0.; } else{ @@ -37,11 +37,10 @@ eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2]; - /*octahedral strain rate*/ - epso = 0.; - for(int i=0;i<3;i++) for(int j=0;j<3;j++) epso += eps[i][j]*eps[i][j]; - if(epso==0.) epso=1.e-14; - epso=sqrt(epso)/sqrt(3.); + /*effective strain rate*/ + epseff = 0.; + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + 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]); - /*Compute the shear strain rate on the non ratating shear plane*/ + /*Compute the shear strain rate on the non rotating shear plane*/ epsprime[0]=0.; epsprime[1]=0.; @@ -74,5 +73,5 @@ /*Assign output pointers*/ - *pepso=epso; + *pepseff=epseff; *pepsprime_norm=epsprime_norm; }/*}}}*/ @@ -114,5 +113,5 @@ omega[0] = 0.; omega[1] = 0.; - omega[2] = 1.; + omega[2] = 0.; } else{ @@ -123,7 +122,14 @@ }/*}}}*/ -IssmDouble EstarLambdaS(IssmDouble epso, IssmDouble epsprime_norm){/*{{{*/ - _assert_(epso>0.); +IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm){/*{{{*/ + IssmDouble lambdas; + _assert_(epsprime_norm>=0.); - return sqrt(2*epsprime_norm*epsprime_norm/(3*epso*epso)); + if(epseff==0.){ + lambdas=0.; + } + else{ + lambdas=sqrt(epsprime_norm*epsprime_norm/(epseff*epseff)); + } + return lambdas; }/*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 21380) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 21381) @@ -15,7 +15,7 @@ IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat); // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); -IssmDouble EstarLambdaS(IssmDouble epso, IssmDouble epsprime_norm); +IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm); void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag); -void EstarStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag); +void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag); IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21380) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21381) @@ -214,6 +214,4 @@ MaterialsRheologyLawEnum, MaterialsRheologyNEnum, - MaterialsRheologyKoEnum, - MaterialsRheologyKobarEnum, MaterialsRheologyEcEnum, MaterialsRheologyEcbarEnum, Index: /issm/trunk-jpl/src/m/classes/matestar.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matestar.m (revision 21380) +++ /issm/trunk-jpl/src/m/classes/matestar.m (revision 21381) @@ -18,5 +18,5 @@ mixed_layer_capacity = 0.; thermal_exchange_velocity = 0.; - rheology_ko = NaN; + rheology_B = NaN; rheology_Ec = NaN; rheology_Es = NaN; @@ -35,5 +35,5 @@ methods function self = extrude(self,md) % {{{ - self.rheology_ko=project3d(md,'vector',self.rheology_ko,'type','node'); + self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node'); self.rheology_Ec=project3d(md,'vector',self.rheology_Ec,'type','node'); self.rheology_Es=project3d(md,'vector',self.rheology_Es,'type','node'); @@ -58,5 +58,5 @@ self.rheology_Es = ones(size(inputstruct.rheology_B)); self.rheology_Ec = ones(size(inputstruct.rheology_B)); - self.rheology_ko = inputstruct.rheology_B.^(-3)*3/2; + self.rheology_B = inputstruct.rheology_B; end otherwise @@ -121,5 +121,5 @@ md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); md = checkfield(md,'fieldname','materials.mu_water','>',0); - md = checkfield(md,'fieldname','materials.rheology_ko','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); md = checkfield(md,'fieldname','materials.rheology_Ec','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); md = checkfield(md,'fieldname','materials.rheology_Es','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); @@ -152,5 +152,5 @@ fielddisplay(self,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); fielddisplay(self,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); - fielddisplay(self,'rheology_ko','octahedral flow law parameter [s^-1 Pa^-2]'); + fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/3)]'); fielddisplay(self,'rheology_Ec','compressive enhancement factor'); fielddisplay(self,'rheology_Es','shear enhancement factor'); @@ -176,5 +176,5 @@ WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double'); WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_ko','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Ec','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Es','format','DoubleMat','mattype',1); @@ -202,5 +202,5 @@ writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity); writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity); - writejs1Darray(fid,[modelname '.materials.rheology_ko'],self.rheology_ko); + writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B); writejs1Darray(fid,[modelname '.materials.rheology_Ec'],self.rheology_Ec); writejs1Darray(fid,[modelname '.materials.rheology_Es'],self.rheology_Es);