Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 13118) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 13119) @@ -100,4 +100,6 @@ MaterialsRheologyLawEnum, MaterialsRheologyNEnum, + MaterialsRheologyZEnum, + MaterialsRheologyZbarEnum, MaterialsRhoIceEnum, MaterialsRhoWaterEnum, @@ -361,4 +363,6 @@ TemperatureOldEnum, TemperaturePicardEnum, + TemperatureSurfaceEnum, + TemperatureBasalEnum, ThicknessAbsMisfitEnum, TypeEnum, Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumToModelField.cpp =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumToModelField.cpp (revision 13119) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumToModelField.cpp (revision 13119) @@ -0,0 +1,30 @@ +/*\file EnumToModelField.cpp: +* \brief: output string associated with enum, that corresponds to a model field +* for example: ThicknessEnum corresponds to model field thickness +* FrictionCoefficientEnum corresponds to model field drag +*/ + +#include "../shared/shared.h" +#include "../include/include.h" +#include "./EnumDefinitions.h" + +const char* EnumToModelField(int en){ + + switch(en){ + + case ThicknessEnum : return "thickness"; + case FrictionCoefficientEnum : return "drag_coefficient"; + case MaterialsRheologyBEnum : return "rheology_B"; + case MaterialsRheologyBbarEnum : return "rheology_B"; + case MaterialsRheologyZEnum : return "rheology_Z"; + case MaterialsRheologyZbarEnum : return "rheology_Z"; + case BalancethicknessThickeningRateEnum: return "dhdt"; + case VxEnum : return "vx"; + case InversionVxObsEnum : return "vx_obs"; + case VyEnum : return "vy"; + case InversionVyObsEnum : return "vy_obs"; + case BasalforcingsMeltingRateEnum : return "basal_melting_rate"; + case SurfaceforcingsMassBalanceEnum : return "surface_mass_balance"; + default : _error_("No model field is associated to enum %s",EnumToStringx(en)); + } +} Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13119) @@ -1555,4 +1555,5 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyBEnum); + else if (enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyZEnum); else input=this->inputs->GetInput(enum_type); //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in penta->inputs"); why error out? if the requested input does not exist, we should still @@ -1667,4 +1668,5 @@ break; case MaterialsRheologyBbarEnum: + case MaterialsRheologyZbarEnum: /*Matice will take care of it*/ break; default: @@ -1962,5 +1964,10 @@ /*update input*/ - this->inputs->AddInput(new PentaP1Input(name,values)); + if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ + matice->inputs->AddInput(new PentaP1Input(name,values)); + } + else{ + this->inputs->AddInput(new PentaP1Input(name,values)); + } return; break; @@ -3230,4 +3237,7 @@ *presponse=this->matice->GetBbar(); break; + case MaterialsRheologyZbarEnum: + *presponse=this->matice->GetZbar(); + break; case VelEnum: { @@ -4460,4 +4470,9 @@ input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); } + else if(enum_type==MaterialsRheologyZbarEnum){ + if(!IsOnBed()) return; + input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); + } + else{ input=inputs->GetInput(enum_type); @@ -4477,4 +4492,7 @@ if(enum_type==MaterialsRheologyBbarEnum){ input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); + } + else if(enum_type==MaterialsRheologyZbarEnum){ + input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); } else{ @@ -4496,4 +4514,7 @@ if(enum_type==MaterialsRheologyBbarEnum){ input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); + } + else if(enum_type==MaterialsRheologyZbarEnum){ + input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); } else{ @@ -4540,4 +4561,7 @@ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + /*Depth Averaging Z*/ + this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + /*Call Tria function*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. @@ -4547,4 +4571,7 @@ /*Delete B averaged*/ this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + + /*Delete Z averaged*/ + this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); /*clean up and return*/ @@ -5111,4 +5138,8 @@ input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input); } + else if(control_type[i]==MaterialsRheologyZbarEnum){ + if (!IsOnBed()) goto cleanup_and_return; + input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input); + } else{ input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input); @@ -5123,4 +5154,7 @@ if(control_type[i]==MaterialsRheologyBbarEnum){ this->InputExtrude(MaterialsRheologyBEnum,MaterialsEnum); + } + else if(control_type[i]==MaterialsRheologyZbarEnum){ + this->InputExtrude(MaterialsRheologyZEnum,MaterialsEnum); } } @@ -6304,4 +6338,7 @@ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + /*Depth Averaging Z*/ + this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + /*Call Tria function*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. @@ -6311,4 +6348,7 @@ /*Delete B averaged*/ this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + + /*Delete Z averaged*/ + this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); /*clean up and return*/ @@ -7860,4 +7900,7 @@ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + /*Depth Averaging Z*/ + this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + /*Call Tria function*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. @@ -7867,4 +7910,7 @@ /*Delete B averaged*/ this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + + /*Delete Z averaged*/ + this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); /*clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13119) @@ -491,4 +491,5 @@ /*Initialize Element matrix*/ ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); @@ -1395,5 +1396,5 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); + if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type); else input=this->inputs->GetInput(enum_type); //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs"); @@ -1501,4 +1502,5 @@ break; case MaterialsRheologyBbarEnum: + case MaterialsRheologyZbarEnum: /*Matice will take care of it*/ break; default: @@ -1689,5 +1691,5 @@ /*update input*/ - if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){ + if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ matice->inputs->AddInput(new TriaP1Input(name,values)); } @@ -1840,4 +1842,5 @@ name==FrictionCoefficientEnum || name==MaterialsRheologyBbarEnum || + name==MaterialsRheologyZbarEnum || name==GradientEnum || name==OldGradientEnum || @@ -2791,4 +2794,7 @@ *presponse=this->matice->GetBbar(); break; + case MaterialsRheologyZbarEnum: + *presponse=this->matice->GetZbar(); + break; case VelEnum:{ @@ -3402,5 +3408,5 @@ for(int i=0;iinputs->GetInput(control_type[i]); _assert_(input); } @@ -3429,5 +3435,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3447,5 +3453,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3466,5 +3472,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3496,4 +3502,7 @@ case MaterialsRheologyBbarEnum: GradjBMacAyeal(gradient,control_index); + break; + case MaterialsRheologyZbarEnum: + GradjZMacAyeal(gradient,control_index); break; case BalancethicknessThickeningRateEnum: @@ -3583,4 +3592,44 @@ } /*}}}*/ +/*FUNCTION Tria::GradjZGradient{{{*/ +void Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){ + + int i,ig; + int doflist1[NUMVERTICES]; + IssmDouble Jdet,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dk[NDOF2]; + IssmDouble grade_g[NUMVERTICES]={0.0}; + GaussTria *gauss=NULL; + + /*Retrieve all inputs we will be needing: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GradientIndexing(&doflist1[0],control_index); + Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + weights_input->GetInputValue(&weight,gauss,weight_index); + + /*Build alpha_complement_list: */ + rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]); + } + gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ /*FUNCTION Tria::GradjBMacAyeal{{{*/ void Tria::GradjBMacAyeal(Vector* gradient,int control_index){ @@ -3624,4 +3673,61 @@ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + /*standard gradient dJ/dki*/ + for (i=0;iweight*basis[i]; + } + + gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL); + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjZMacAyeal{{{*/ +void Tria::GradjZMacAyeal(Vector* gradient,int control_index){ + + /*Intermediaries*/ + int i,ig; + int doflist[NUMVERTICES]; + IssmDouble vx,vy,lambda,mu,thickness,Jdet; + IssmDouble viscosity_complement; + IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3],epsilon[3]; + IssmDouble grad[NUMVERTICES]={0.0}; + GaussTria *gauss = NULL; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GradientIndexing(&doflist[0],control_index); + + /*Retrieve all inputs*/ + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(4); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness,gauss); + rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss); + vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); + adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); + + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); @@ -5109,5 +5215,5 @@ /*Get input (either in element or material)*/ - if(control_enum==MaterialsRheologyBbarEnum){ + if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); } @@ -5144,5 +5250,5 @@ new_input = new TriaP1Input(control_enum,values); - if(control_enum==MaterialsRheologyBbarEnum){ + if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); } Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 13118) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 13119) @@ -143,5 +143,7 @@ void Gradj(Vector* gradient,int control_type,int control_index); void GradjBGradient(Vector* gradient,int weight_index,int control_index); + void GradjZGradient(Vector* gradient,int weight_index,int control_index); void GradjBMacAyeal(Vector* gradient,int control_index); + void GradjZMacAyeal(Vector* gradient,int control_index); void GradjDragMacAyeal(Vector* gradient,int control_index); void GradjDragStokes(Vector* gradient,int control_index); Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp (revision 13119) @@ -172,4 +172,24 @@ } /*}}}*/ +/*FUNCTION Matice::GetZ {{{*/ +IssmDouble Matice::GetZ(){ + + /*Output*/ + IssmDouble Z; + + inputs->GetInputAverage(&Z,MaterialsRheologyZEnum); + return Z; +} +/*}}}*/ +/*FUNCTION Matice::GetZbar {{{*/ +IssmDouble Matice::GetZbar(){ + + /*Output*/ + IssmDouble Zbar; + + inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum); + return Zbar; +} +/*}}}*/ /*FUNCTION Matice::GetVectorFromInputs{{{*/ void Matice::GetVectorFromInputs(Vector* vector,int input_enum){ @@ -207,5 +227,5 @@ void Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ /*From a string tensor and a material object, return viscosity, using Glen's flow law. - B + Z * B viscosity= ------------------------------------------------------------------- 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -226,9 +246,11 @@ /*Intermediary: */ IssmDouble A,e; - IssmDouble B,n; + IssmDouble Btmp,B,n,Z; /*Get B and n*/ - B=GetBbar(); + Btmp=GetBbar(); + Z=GetZbar(); n=GetN(); + B=Z*Btmp; if (n==1){ @@ -292,9 +314,10 @@ /*Intermediaries: */ IssmDouble A,e; - IssmDouble B,n; - - /*Get B and n*/ - B=GetB(); + IssmDouble B,n,Z; + + /*Get B, Z and n*/ n=GetN(); + Z=GetZ(); + B=Z*GetB(); if (n==1){ @@ -362,11 +385,12 @@ /*Intermediaries: */ IssmDouble A,e; - IssmDouble B,n; + IssmDouble B,n,Z; IssmDouble eps0; /*Get B and n*/ eps0=pow((IssmDouble)10,(IssmDouble)-27); - B=GetB(); n=GetN(); + Z=GetZ(); + B=Z*GetB(); if (n==1){ @@ -452,4 +476,60 @@ viscosity_complement=1/(2*pow(A,e)); + } + } + else{ + viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); + } + + /*Checks in debugging mode*/ + _assert_(B>0); + _assert_(n>0); + _assert_(viscosity_complement>0); + + /*Return: */ + *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosityZComplement {{{*/ +void Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ + /*Return viscosity derivative for control method d(mu)/dZ: + * + * B + * dviscosity= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * If epsilon is NULL, it means this is the first time Gradjb is being run, and we + * return mu20, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity_complement; + + /*input strain rate: */ + IssmDouble exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + IssmDouble A,e; + IssmDouble B,n; + + /*Get B and n*/ + B=GetBbar(); + n=GetN(); + + if(epsilon){ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: mu2=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity_complement for 0 shear areas: */ + viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); + } + else{ + e=(n-1)/(2*n); + + viscosity_complement=B/(2*pow(A,e)); } } @@ -701,4 +781,10 @@ } + /*Get Z*/ + if (iomodel->Data(MaterialsRheologyZEnum)) { + for(i=0;iData(MaterialsRheologyZEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); + } + /*Control Inputs*/ #ifdef _HAVE_CONTROL_ @@ -715,4 +801,14 @@ } break; + case MaterialsRheologyZbarEnum: + if (iomodel->Data(MaterialsRheologyZEnum)){ + _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + } } @@ -741,4 +837,10 @@ for(i=0;iData(MaterialsRheologyNEnum)[index]; this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); + } + + /*Get Z*/ + if (iomodel->Data(MaterialsRheologyZEnum)) { + for(i=0;iData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs)); } @@ -757,4 +859,13 @@ } break; + case MaterialsRheologyZbarEnum: + if (iomodel->Data(MaterialsRheologyZEnum)){ + _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; } } @@ -775,5 +886,7 @@ name==MaterialsRheologyBEnum || name==MaterialsRheologyBbarEnum || - name==MaterialsRheologyNEnum + name==MaterialsRheologyNEnum || + name==MaterialsRheologyZEnum || + name==MaterialsRheologyZbarEnum ){ return true; Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h (revision 13118) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h (revision 13119) @@ -58,16 +58,19 @@ /*}}}*/ /*Matice Numerics: {{{*/ - void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); - void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); - void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); - void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); - void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); + void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); + void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); + void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); + void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); IssmDouble GetA(); IssmDouble GetB(); IssmDouble GetBbar(); IssmDouble GetN(); - bool IsInput(int name); + IssmDouble GetZ(); + IssmDouble GetZbar(); + bool IsInput(int name); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 13119) @@ -105,4 +105,6 @@ case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; case MaterialsRheologyNEnum : return "MaterialsRheologyN"; + case MaterialsRheologyZEnum : return "MaterialsRheologyZ"; + case MaterialsRheologyZbarEnum : return "MaterialsRheologyZbar"; case MaterialsRhoIceEnum : return "MaterialsRhoIce"; case MaterialsRhoWaterEnum : return "MaterialsRhoWater"; @@ -352,4 +354,6 @@ case TemperatureOldEnum : return "TemperatureOld"; case TemperaturePicardEnum : return "TemperaturePicard"; + case TemperatureSurfaceEnum : return "TemperatureSurface"; + case TemperatureBasalEnum : return "TemperatureBasal"; case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit"; case TypeEnum : return "Type"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 13119) @@ -48,4 +48,5 @@ case FrictionCoefficientEnum: iomodel->FetchData(1,FrictionCoefficientEnum); break; case MaterialsRheologyBbarEnum: iomodel->FetchData(1,MaterialsRheologyBEnum); break; + case MaterialsRheologyZbarEnum: iomodel->FetchData(1,MaterialsRheologyZEnum); break; default: _error_("Control " << EnumToStringx(reCast(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); } @@ -66,4 +67,4 @@ /*Free data: */ - iomodel->DeleteData(1+4+5,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum); + iomodel->DeleteData(1+4+6,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum,MaterialsRheologyZEnum); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 13119) @@ -45,5 +45,5 @@ /*Fetch data needed: */ - iomodel->FetchData(4,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum); + iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum); #ifdef _HAVE_3D_ if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum); @@ -67,8 +67,9 @@ /*Free data: */ - iomodel->DeleteData(9,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum, - MaterialsRheologyBEnum,MaterialsRheologyNEnum,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum); + iomodel->DeleteData(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum, + MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum,InversionControlParametersEnum,InversionMinParametersEnum, + InversionMaxParametersEnum); - /*Add new constrant material property tgo materials, at the end: */ + /*Add new constant material property to materials, at the end: */ materials->AddObject(new Matpar(numberofelements+1,iomodel));//put it at the end of the materials @@ -81,5 +82,5 @@ for (i=0;imy_vertices[i]){ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 13119) @@ -61,4 +61,5 @@ iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum); iomodel->FetchDataToInput(elements,VxEnum); iomodel->FetchDataToInput(elements,VyEnum); Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 13119) @@ -47,5 +47,5 @@ } - /*Assign output pointers: */ + /*Assign output pointers:*/ *puf=uf; } Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 13118) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 13119) @@ -106,4 +106,6 @@ else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; + else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum; + else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum; else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum; else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum; @@ -136,10 +138,10 @@ else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum; else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum; - else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum; - else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum; + if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum; + else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum; + else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum; else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum; else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum; @@ -259,10 +261,10 @@ else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else if (strcmp(name,"Contour")==0) return ContourEnum; - else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; - else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum; + if (strcmp(name,"ControlInput")==0) return ControlInputEnum; + else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; + else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; @@ -359,4 +361,6 @@ else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum; else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; + else if (strcmp(name,"TemperatureSurface")==0) return TemperatureSurfaceEnum; + else if (strcmp(name,"TemperatureBasal")==0) return TemperatureBasalEnum; else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; else if (strcmp(name,"Type")==0) return TypeEnum; @@ -380,12 +384,12 @@ else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; - else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum; else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; + else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 13118) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 13119) @@ -12,4 +12,6 @@ pressure = NaN; temperature = NaN; + surfacetemp = NaN; + basaltemp = NaN; watercolumn = NaN; waterfraction = NaN; @@ -68,4 +70,6 @@ fielddisplay(obj,'pressure','pressure field'); fielddisplay(obj,'temperature','temperature in Kelvins'); + fielddisplay(obj,'surfacetemp','surface temperature in Kelvins'); + fielddisplay(obj,'basaltemp','basal temperature in Kelvins'); fielddisplay(obj,'watercolumn','thickness of subglacial water'); fielddisplay(obj,'waterfraction','fraction of water in the ice'); @@ -73,11 +77,13 @@ end % }}} function marshall(obj,fid) % {{{ - WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum()); - WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum()); - WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum()); - WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum()); - WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum()); - WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum()); - WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum()); + WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum); + WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum); + WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum); + WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum); + WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum); + WriteData(fid,'data',obj.surfacetemp,'format','DoubleMat','mattype',1,'enum',TemperatureSurfaceEnum); + WriteData(fid,'data',obj.basaltemp,'format','DoubleMat','mattype',1,'enum',TemperatureBasalEnum); + WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum); + WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum); end % }}} end Index: /issm/trunk-jpl/src/m/classes/inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.m (revision 13118) +++ /issm/trunk-jpl/src/m/classes/inversion.m (revision 13119) @@ -83,17 +83,17 @@ num_costfunc=size(md.inversion.cost_functions,2); - md = checkfield(md,'inversion.iscontrol','values',[0 1]); - md = checkfield(md,'inversion.tao','values',[0 1]); - md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]); - md = checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'}); - md = checkfield(md,'inversion.nsteps','numel',[1],'>=',1); - md = checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0); - md = checkfield(md,'inversion.step_threshold','size',[md.inversion.nsteps 1]); - md = checkfield(md,'inversion.cost_functions','size',[md.inversion.nsteps num_costfunc],'values',[101:105 201 501:505]); - md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0); - md = checkfield(md,'inversion.gradient_only','values',[0 1]); - md = checkfield(md,'inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]); - md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]); - md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]); + checkfield(md,'inversion.iscontrol','values',[0 1]); + checkfield(md,'inversion.tao','values',[0 1]); + checkfield(md,'inversion.incomplete_adjoint','values',[0 1]); + checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy'}); + checkfield(md,'inversion.nsteps','numel',1,'>=',1); + checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0); + checkfield(md,'inversion.step_threshold','size',[md.inversion.nsteps 1]); + checkfield(md,'inversion.cost_functions','size',[md.inversion.nsteps num_costfunc],'values',[101:105 201 501:503]); + checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0); + checkfield(md,'inversion.gradient_only','values',[0 1]); + checkfield(md,'inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]); + checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]); + checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]); if solution==BalancethicknessSolutionEnum() Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 13118) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 13119) @@ -19,4 +19,5 @@ rheology_B = NaN; rheology_n = NaN; + rheology_Z = NaN; rheology_law = ''; end @@ -76,4 +77,5 @@ md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'}); end % }}} @@ -94,4 +96,5 @@ fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]'); fielddisplay(obj,'rheology_n','Glen''s flow law exponent'); + fielddisplay(obj,'rheology_Z','rheology multiplier'); fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius'''); end % }}} @@ -110,4 +113,5 @@ WriteData(fid,'object',obj,'fieldname','rheology_B','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','rheology_n','format','DoubleMat','mattype',2); + WriteData(fid,'object',obj,'fieldname','rheology_Z','format','DoubleMat','mattype',1); WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer'); end % }}} Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 13118) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 13119) @@ -168,4 +168,5 @@ md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); + md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); %special for thermal modeling: @@ -716,4 +717,5 @@ md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); + md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); %parameters @@ -829,6 +831,8 @@ if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end + if isfield(structmd,'Z'), md.materials.rheology_Z=structmd.Z; end if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end + if isfield(structmd,'rheology_Z'), md.materials.rheology_Z=structmd.rheology_Z; end if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end Index: /issm/trunk-jpl/src/m/contrib/hack/tres.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/hack/tres.m (revision 13118) +++ /issm/trunk-jpl/src/m/contrib/hack/tres.m (revision 13119) @@ -105,4 +105,6 @@ case MaterialsRheologyBEnum(), string='rheology_B'; return case MaterialsRheologyBbarEnum(), string='rheology_B'; return + case MaterialsRheologyZEnum(), string='rheology_Z'; return + case MaterialsRheologyZbarEnum(), string='rheology_Z'; return case BalancethicknessThickeningRateEnum(), string='dhdt'; return case VxEnum(), string='vx'; return Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 13118) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 13119) @@ -889,4 +889,24 @@ return StringToEnum('MaterialsRheologyN')[0] +def MaterialsRheologyZEnum(): + """ + MATERIALSRHEOLOGYZENUM - Enum of MaterialsRheologyZ + + Usage: + macro=MaterialsRheologyZEnum() + """ + + return StringToEnum('MaterialsRheologyZ')[0] + +def MaterialsRheologyZbarEnum(): + """ + MATERIALSRHEOLOGYZBARENUM - Enum of MaterialsRheologyZbar + + Usage: + macro=MaterialsRheologyZbarEnum() + """ + + return StringToEnum('MaterialsRheologyZbar')[0] + def MaterialsRhoIceEnum(): """ @@ -3359,4 +3379,24 @@ return StringToEnum('TemperaturePicard')[0] +def TemperatureSurfaceEnum(): + """ + TEMPERATURESURFACEENUM - Enum of TemperatureSurface + + Usage: + macro=TemperatureSurfaceEnum() + """ + + return StringToEnum('TemperatureSurface')[0] + +def TemperatureBasalEnum(): + """ + TEMPERATUREBASALENUM - Enum of TemperatureBasal + + Usage: + macro=TemperatureBasalEnum() + """ + + return StringToEnum('TemperatureBasal')[0] + def ThicknessAbsMisfitEnum(): """ @@ -4607,4 +4647,4 @@ """ - return 459 - + return 463 + Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 13118) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 13119) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=459; +macro=463; Index: /issm/trunk-jpl/src/m/enum/TemperatureBasalEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/TemperatureBasalEnum.m (revision 13119) +++ /issm/trunk-jpl/src/m/enum/TemperatureBasalEnum.m (revision 13119) @@ -0,0 +1,11 @@ +function macro=TemperatureBasalEnum() +%TEMPERATUREBASALENUM - Enum of TemperatureBasal +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=TemperatureBasalEnum() + +macro=StringToEnum('TemperatureBasal'); Index: /issm/trunk-jpl/src/m/enum/TemperatureSurfaceEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/TemperatureSurfaceEnum.m (revision 13119) +++ /issm/trunk-jpl/src/m/enum/TemperatureSurfaceEnum.m (revision 13119) @@ -0,0 +1,11 @@ +function macro=TemperatureSurfaceEnum() +%TEMPERATURESURFACEENUM - Enum of TemperatureSurface +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=TemperatureSurfaceEnum() + +macro=StringToEnum('TemperatureSurface'); Index: /issm/trunk-jpl/test/NightlyRun/IdToName.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/IdToName.m (revision 13118) +++ /issm/trunk-jpl/test/NightlyRun/IdToName.m (revision 13119) @@ -65,4 +65,7 @@ case 236, name='SquareShelfTranIspddIsdeltaM2d'; case 237, name='SquareShelfTranIspddIsdeltaM3d'; + case 270, name='SquareShelfDiagM2dDamage'; + case 272, name='SquareShelfCMZM2dDamage'; + case 274, name='SquareShelfDiagM2dDamageRift'; case 301, name='SquareSheetConstrainedDiagM2d'; case 302, name='SquareSheetConstrainedDiagH2d'; Index: /issm/trunk-jpl/test/NightlyRun/README =================================================================== --- /issm/trunk-jpl/test/NightlyRun/README (revision 13118) +++ /issm/trunk-jpl/test/NightlyRun/README (revision 13119) @@ -31,5 +31,5 @@ 6: 79North Add the id and testname in IdToName.m (incresing order) -We try not to create to many .par and .exp files, so try to use the existing ones as much as possible. +We try not to create too many .par and .exp files, so try to use the existing ones as much as possible. To modify some characteristics, do it in the testxxx.m file. The testxxx_nightly.m is used to define the parameters you want to check in the nightlyruns. Index: /issm/trunk-jpl/test/NightlyRun/test270.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test270.m (revision 13119) +++ /issm/trunk-jpl/test/NightlyRun/test270.m (revision 13119) @@ -0,0 +1,17 @@ +md=triangle(model,'../Exp/Square.exp',150000); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md.materials.rheology_Z=0.5*ones(md.mesh.numberofvertices,1); +md=setflowequation(md,'macayeal','all'); +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,DiagnosticSolutionEnum); + +%Fields and tolerances to track changes +field_names ={'Vx','Vy','Vel','Pressure'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.DiagnosticSolution.Vx),... + (md.results.DiagnosticSolution.Vy),... + (md.results.DiagnosticSolution.Vel),... + (md.results.DiagnosticSolution.Pressure),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test272.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test272.m (revision 13119) +++ /issm/trunk-jpl/test/NightlyRun/test272.m (revision 13119) @@ -0,0 +1,35 @@ +md=triangle(model,'../Exp/Square.exp',150000); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md.materials.rheology_Z=0.5*ones(md.mesh.numberofvertices,1); +md=setflowequation(md,'macayeal','all'); + +%control parameters +md.inversion.iscontrol=1; +md.inversion.control_parameters={'MaterialsRheologyZbar'}; +md.inversion.min_parameters=10^-13*ones(md.mesh.numberofvertices,1); +md.inversion.max_parameters=ones(md.mesh.numberofvertices,1); +md.inversion.nsteps=2; +md.inversion.cost_functions=101*ones(md.inversion.nsteps,1); +md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); +md.inversion.gradient_scaling=0.9*ones(md.inversion.nsteps,1); +md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1); +md.inversion.step_threshold=0.99*ones(md.inversion.nsteps,1); +md.inversion.vx_obs=md.initialization.vx; +md.inversion.vy_obs=md.initialization.vy; + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,DiagnosticSolutionEnum); + +%Fields and tolerances to track changes +field_names ={'Gradient' 'Misfits' 'MaterialsRheologyZbar' 'Pressure' 'Vel' 'Vx' 'Vy'}; +field_tolerances={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12}; +field_values={... + (md.results.DiagnosticSolution.Gradient1),... + md.results.DiagnosticSolution.J,... + (md.results.DiagnosticSolution.MaterialsRheologyZbar),... + (md.results.DiagnosticSolution.Pressure),... + (md.results.DiagnosticSolution.Vel),... + (md.results.DiagnosticSolution.Vx),... + (md.results.DiagnosticSolution.Vy) +}; Index: /issm/trunk-jpl/test/NightlyRun/test274.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test274.m (revision 13119) +++ /issm/trunk-jpl/test/NightlyRun/test274.m (revision 13119) @@ -0,0 +1,19 @@ +md=triangle(model,'../Exp/SquareHole.exp','../Exp/Rifts.exp',50000); +md=meshprocessrifts(md,'../Exp/Square.exp'); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md.materials.rheology_Z=0.5*ones(md.mesh.numberofvertices,1); +md=setflowequation(md,'macayeal','all'); + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,DiagnosticSolutionEnum); + +%Fields and tolerances to track changes +field_names ={'Vx','Vy','Vel','Pressure'}; +field_tolerances={1e-11,1e-11,1e-11,1e-11}; +field_values={... + (md.results.DiagnosticSolution.Vx),... + (md.results.DiagnosticSolution.Vy),... + (md.results.DiagnosticSolution.Vel),... + (md.results.DiagnosticSolution.Pressure),... + }; Index: /issm/trunk-jpl/test/Par/79North.par =================================================================== --- /issm/trunk-jpl/test/Par/79North.par (revision 13118) +++ /issm/trunk-jpl/test/Par/79North.par (revision 13119) @@ -13,4 +13,5 @@ md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); Index: /issm/trunk-jpl/test/Par/ISMIPA.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPA.par (revision 13118) +++ /issm/trunk-jpl/test/Par/ISMIPA.par (revision 13119) @@ -17,4 +17,5 @@ md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' boundary conditions for diagnostic model'); Index: /issm/trunk-jpl/test/Par/ISMIPB.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPB.par (revision 13118) +++ /issm/trunk-jpl/test/Par/ISMIPB.par (revision 13119) @@ -17,4 +17,5 @@ md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' boundary conditions for diagnostic model'); Index: /issm/trunk-jpl/test/Par/ISMIPC.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPC.par (revision 13118) +++ /issm/trunk-jpl/test/Par/ISMIPC.par (revision 13119) @@ -18,4 +18,5 @@ md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' boundary conditions for diagnostic model: '); Index: /issm/trunk-jpl/test/Par/ISMIPD.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPD.par (revision 13118) +++ /issm/trunk-jpl/test/Par/ISMIPD.par (revision 13119) @@ -17,4 +17,5 @@ md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' boundary conditions for diagnostic model: '); Index: /issm/trunk-jpl/test/Par/ISMIPE.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPE.par (revision 13118) +++ /issm/trunk-jpl/test/Par/ISMIPE.par (revision 13119) @@ -26,4 +26,5 @@ md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' boundary conditions for diagnostic model: '); Index: /issm/trunk-jpl/test/Par/ISMIPF.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPF.par (revision 13118) +++ /issm/trunk-jpl/test/Par/ISMIPF.par (revision 13119) @@ -16,4 +16,5 @@ md.materials.rheology_B=1.4734*10^14*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=1*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); md.materials.rheology_law='None'; Index: /issm/trunk-jpl/test/Par/Pig.par =================================================================== --- /issm/trunk-jpl/test/Par/Pig.par (revision 13118) +++ /issm/trunk-jpl/test/Par/Pig.par (revision 13119) @@ -17,4 +17,5 @@ md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); md.initialization.temperature=md.initialization.temperature; Index: /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par (revision 13118) +++ /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par (revision 13119) @@ -20,4 +20,5 @@ md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' creating surface mass balance'); Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 13118) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 13119) @@ -60,4 +60,5 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); %Surface mass balance and basal melting Index: /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par (revision 13118) +++ /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par (revision 13119) @@ -25,4 +25,5 @@ md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' creating surface mass balance'); Index: /issm/trunk-jpl/test/Par/SquareEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareEISMINT.par (revision 13118) +++ /issm/trunk-jpl/test/Par/SquareEISMINT.par (revision 13119) @@ -27,4 +27,5 @@ md.materials.rheology_B=1.7687*10^8*ones(md.mesh.numberofvertices,1); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' creating surface mass balance'); Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.par (revision 13118) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.par (revision 13119) @@ -22,4 +22,5 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); %Friction Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetShelf.par (revision 13118) +++ /issm/trunk-jpl/test/Par/SquareSheetShelf.par (revision 13119) @@ -25,4 +25,5 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); %Accumulation and melting Index: /issm/trunk-jpl/test/Par/SquareShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf.par (revision 13118) +++ /issm/trunk-jpl/test/Par/SquareShelf.par (revision 13119) @@ -22,4 +22,5 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); %Friction Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelfConstrained.par (revision 13118) +++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.par (revision 13119) @@ -22,4 +22,5 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); %Surface mass balance and basal melting Index: /issm/trunk-jpl/test/Par/SquareThermal.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareThermal.par (revision 13118) +++ /issm/trunk-jpl/test/Par/SquareThermal.par (revision 13119) @@ -28,4 +28,5 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_Z=ones(md.mesh.numberofvertices,1); disp(' creating surface mass balance');