Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16166) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16167) @@ -100,6 +100,4 @@ ./classes/Materials/Matice.h\ ./classes/Materials/Matice.cpp\ - ./classes/Materials/Matdamageice.h\ - ./classes/Materials/Matdamageice.cpp\ ./classes/Materials/Matpar.h\ ./classes/Materials/Matpar.cpp\ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16167) @@ -2142,10 +2142,7 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==MaterialsRheologyBbarEnum) - input=this->material->inputs->GetInput(MaterialsRheologyBEnum); - else if (enum_type==MaterialsRheologyZbarEnum) - input=this->material->inputs->GetInput(MaterialsRheologyZEnum); - else - input=this->inputs->GetInput(enum_type); + if (enum_type==MaterialsRheologyBbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyBEnum); + else if (enum_type==DamageDbarEnum) input=this->material->inputs->GetInput(DamageDEnum); + 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 //try and output whatever we can instead of just failing. @@ -2258,7 +2255,7 @@ } break; - case MaterialsRheologyBbarEnum: - case MaterialsRheologyZbarEnum: - /*Material will take care of it*/ break; + /*Material will take care of it*/ + case MaterialsRheologyBbarEnum: break; + case DamageDbarEnum:break; default: _error_("Control " << EnumToStringx(reCast(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); @@ -2667,10 +2664,5 @@ } /*update input*/ - if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new PentaInput(name,values,P1Enum)); - } - else{ - this->inputs->AddInput(new PentaInput(name,values,P1Enum)); - } + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); return; @@ -2680,10 +2672,5 @@ } /*update input*/ - if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new PentaInput(name,values,P1Enum)); - } - else{ - this->inputs->AddInput(new PentaInput(name,values,P1Enum)); - } + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); return; @@ -3980,6 +3967,6 @@ *presponse=this->material->GetBbar(); break; - case MaterialsRheologyZbarEnum: - *presponse=this->material->GetZbar(); + case DamageDbarEnum: + *presponse=this->material->GetDbar(); break; case VelEnum: @@ -5433,9 +5420,8 @@ input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); } - else if(enum_type==MaterialsRheologyZbarEnum){ + else if(enum_type==DamageDbarEnum){ if(!IsOnBed()) return; - input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); - } - + input=(Input*)material->inputs->GetInput(DamageDEnum); + } else{ input=inputs->GetInput(enum_type); @@ -5456,6 +5442,6 @@ input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); } - else if(enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); + else if(enum_type==DamageDbarEnum){ + input=(Input*)material->inputs->GetInput(DamageDEnum); } else{ @@ -5478,6 +5464,6 @@ input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); } - else if(enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); + else if(enum_type==DamageDbarEnum){ + input=(Input*)material->inputs->GetInput(DamageDEnum); } else{ @@ -5525,8 +5511,5 @@ case MaticeEnum: this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - break; - case MatdamageiceEnum: - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); break; default: @@ -5541,5 +5524,5 @@ /*Delete averaged fields*/ this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); /*clean up and return*/ @@ -6137,5 +6120,5 @@ /*If on water, skip grad (=0): */ if(NoIceInElement())return; - + /*First deal with ∂/∂alpha(KU-F)*/ switch(control_type){ @@ -6404,4 +6387,5 @@ /*Depth Average B*/ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); /*Collapse element to the base*/ @@ -6412,4 +6396,5 @@ /*delete Average B*/ this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); } /*}}}*/ @@ -6420,6 +6405,7 @@ if (!IsOnBed()) return; - /*Depth Average B*/ + /*Depth Average B and D*/ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); /*Collapse element to the base*/ @@ -6430,4 +6416,5 @@ /*delete Average B*/ this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); } /*}}}*/ /*FUNCTION Penta::GradjBbarFS {{{*/ @@ -6437,6 +6424,7 @@ if (!IsOnBed()) return; - /*Depth Average B*/ + /*Depth Average B and D*/ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); /*Collapse element to the base*/ @@ -6447,4 +6435,5 @@ /*delete Average B*/ this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); } /*}}}*/ /*FUNCTION Penta::InputControlUpdate{{{*/ @@ -6466,7 +6455,7 @@ input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input); } - else if(control_type[i]==MaterialsRheologyZbarEnum){ + else if(control_type[i]==DamageDbarEnum){ if (!IsOnBed()) goto cleanup_and_return; - input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input); + input=(Input*)material->inputs->GetInput(DamageDEnum); _assert_(input); } else{ @@ -6483,6 +6472,6 @@ this->InputExtrude(MaterialsRheologyBEnum,MaterialsEnum); } - else if(control_type[i]==MaterialsRheologyZbarEnum){ - this->InputExtrude(MaterialsRheologyZEnum,MaterialsEnum); + else if(control_type[i]==DamageDbarEnum){ + this->InputExtrude(DamageDEnum,MaterialsEnum); } } @@ -7695,8 +7684,5 @@ case MaticeEnum: this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - break; - case MatdamageiceEnum: - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); break; default: @@ -7711,5 +7697,5 @@ /*Delete averaged fields*/ this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); /*clean up and return*/ @@ -9737,8 +9723,5 @@ case MaticeEnum: this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - break; - case MatdamageiceEnum: - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); break; default: @@ -9753,5 +9736,5 @@ /*Delete averaged inputs*/ this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); /*clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16167) @@ -1561,7 +1561,7 @@ } break; - case MaterialsRheologyBbarEnum: - case DamageDbarEnum: - /*Material will take care of it*/ break; + /*Material will take care of it*/ + case MaterialsRheologyBbarEnum:break; + case DamageDbarEnum: break; default: _error_("Control " << EnumToStringx(reCast(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 16167) @@ -61,4 +61,5 @@ /*FUNCTION ControlInput::Echo {{{*/ void ControlInput::Echo(void){ + printf("I'm here\n"); this->DeepEcho(); } Index: sm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 16166) +++ (revision ) @@ -1,912 +1,0 @@ -/*!\file Matdamageice.c - * \brief: implementation of the Matdamageice object - */ - -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "../classes.h" -#include "shared/shared.h" - -/*Matdamageice constructors and destructor*/ -/*FUNCTION Matdamageice::Matdamageice(){{{*/ -Matdamageice::Matdamageice(){ - this->inputs=NULL; - this->helement=NULL; - return; -} -/*}}}*/ -/*FUNCTION Matdamageice::Matdamageice(int id, int index, IoModel* iomodel, int num_vertices){{{*/ -Matdamageice::Matdamageice(int matice_mid,int index, IoModel* iomodel){ - - /*Intermediaries:*/ - int matice_eid; - - /*Initialize id*/ - this->mid=matice_mid; - - /*Initialize inputs*/ - this->inputs=new Inputs(); - - /*Initialize inputs from IoModel*/ - this->InputUpdateFromIoModel(index,iomodel); - - /*Hooks: */ - matice_eid=index+1; - this->helement=new Hook(&matice_eid,1); - - return; - -} -/*}}}*/ -/*FUNCTION Matdamageice::~Matdamageice(){{{*/ -Matdamageice::~Matdamageice(){ - delete helement; - delete inputs; - return; -} -/*}}}*/ - -/*Object virtual functions definitions:*/ -/*FUNCTION Matdamageice::Echo {{{*/ -void Matdamageice::Echo(void){ - - _printf_("Matdamageice:\n"); - _printf_(" mid: " << mid << "\n"); - _printf_(" inputs:\n"); - inputs->Echo(); - _printf_(" element:\n"); - helement->Echo(); -} -/*}}}*/ -/*FUNCTION Matdamageice::DeepEcho {{{*/ -void Matdamageice::DeepEcho(void){ - - _printf_("Matdamageice:\n"); - _printf_(" mid: " << mid << "\n"); - _printf_(" inputs:\n"); - inputs->DeepEcho(); - _printf_(" element:\n"); - helement->Echo(); -} -/*}}}*/ -/*FUNCTION Matdamageice::Id {{{*/ -int Matdamageice::Id(void){ return mid; } -/*}}}*/ -/*FUNCTION Matdamageice::ObjectEnum{{{*/ -int Matdamageice::ObjectEnum(void){ - - return MatdamageiceEnum; - -} -/*}}}*/ -/*FUNCTION Matdamageice::copy {{{*/ -Object* Matdamageice::copy() { - - /*Output*/ - Matdamageice* matice=NULL; - - /*Initialize output*/ - matice=new Matdamageice(); - - /*copy fields: */ - matice->mid=this->mid; - matice->helement=(Hook*)this->helement->copy(); - if(this->inputs) matice->inputs=(Inputs*)this->inputs->Copy(); - else matice->inputs=new Inputs(); - - return matice; -} -/*}}}*/ - -/*Matdamageice management*/ -/*FUNCTION Matdamageice::Configure {{{*/ -void Matdamageice::Configure(Elements* elementsin){ - - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective - * datasets, using internal ids and offsets hidden in hooks: */ - helement->configure((DataSet*)elementsin); -} -/*}}}*/ -/*FUNCTION Matdamageice::SetCurrentConfiguration {{{*/ -void Matdamageice::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){ - -} -/*}}}*/ -/*FUNCTION Matdamageice::GetA {{{*/ -IssmDouble Matdamageice::GetA(){ - /* - * A = 1/B^n - */ - - IssmDouble B,n; - - inputs->GetInputAverage(&B,MaterialsRheologyBEnum); - n=this->GetN(); - - return pow(B,-n); -} -/*}}}*/ -/*FUNCTION Matdamageice::GetB {{{*/ -IssmDouble Matdamageice::GetB(){ - - /*Output*/ - IssmDouble B; - - inputs->GetInputAverage(&B,MaterialsRheologyBEnum); - return B; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetBbar {{{*/ -IssmDouble Matdamageice::GetBbar(){ - - /*Output*/ - IssmDouble Bbar; - - inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum); - return Bbar; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetN {{{*/ -IssmDouble Matdamageice::GetN(){ - - /*Output*/ - IssmDouble n; - - inputs->GetInputAverage(&n,MaterialsRheologyNEnum); - return n; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetD {{{*/ -IssmDouble Matdamageice::GetD(){ - - /*Output*/ - IssmDouble D; - - inputs->GetInputAverage(&D,DamageDEnum); - return D; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetZ {{{*/ -IssmDouble Matdamageice::GetZ(){ - - /*Output*/ - IssmDouble D; - - inputs->GetInputAverage(&D,DamageDEnum); - return 1/(1-D); -} -/*}}}*/ -/*FUNCTION Matdamageice::GetZbar {{{*/ -IssmDouble Matdamageice::GetZbar(){ - - /*Output*/ - IssmDouble Dbar; - inputs->GetInputAverage(&Dbar,DamageDbarEnum); - return 1/(1-Dbar); -} -/*}}}*/ -/*FUNCTION Matdamageice::GetDbar {{{*/ -IssmDouble Matdamageice::GetDbar(){ - - /*Output*/ - IssmDouble Dbar; - inputs->GetInputAverage(&Dbar,DamageDbarEnum); - return Dbar; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetVectorFromInputs{{{*/ -void Matdamageice::GetVectorFromInputs(Vector* vector,int input_enum){ - - /*Intermediaries*/ - Element *element= NULL; - - /*Recover element*/ - element=(Element*)helement->delivers(); - - /*Check that input_enum is a material input*/ - if (!IsInput(input_enum)) return; - - switch(element->ObjectEnum()){ - - case TriaEnum:{ - - /*Prepare index list*/ - int vertexpidlist[3]; - ((Tria*)element)->GetVertexPidList(&vertexpidlist[0]); - - /*Get input (either in element or material)*/ - Input* input=inputs->GetInput(input_enum); - if(!input) _error_("Input " << EnumToStringx(input_enum) << " not found in material"); - - /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ - input->GetVectorFromInputs(vector,&vertexpidlist[0]);} - break; - - default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet"); - } -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosity2d {{{*/ -void Matdamageice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ - /*From a string tensor and a material object, return viscosity, using Glen's flow law. - (1-D) * B - viscosity= ------------------------------------------------------------------- - 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] - - where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity - vector, and n the flow law exponent. - - If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we - return 10^14, initial viscosity. - */ - - /*output: */ - IssmDouble viscosity; - - /*input strain rate: */ - IssmDouble exx,eyy,exy; - - /*Intermediary: */ - IssmDouble A,e; - IssmDouble Btmp,B,n,D; - - /*Get B and n*/ - Btmp=GetBbar(); - D=GetDbar(); - n=GetN(); - B=(1-D)*Btmp; - - if (n==1){ - /*Viscous behaviour! viscosity=B: */ - viscosity=B/2; - } - else{ - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ - viscosity=0.5*pow(10.,14); - } - else{ - /*Retrive strain rate components: */ - exx=*(epsilon+0); - eyy=*(epsilon+1); - exy=*(epsilon+2); - - /*Build viscosity: viscosity=B/(2*A^e) */ - A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; - if(A==0){ - /*Maxiviscositym viscosity for 0 shear areas: */ - viscosity=2.5*pow(10.,17); - } - else{ - e=(n-1)/(2*n); - viscosity=B/(2*pow(A,e)); - } - } - } - - /*Checks in debugging mode*/ - if(viscosity<=0) _error_("Negative viscosity"); - _assert_(B>0); - _assert_(n>0); - - /*Return: */ - *pviscosity=viscosity; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosity3d {{{*/ -void Matdamageice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){ - - /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: - * - * (1-D)* B - * viscosity3d= ------------------------------------------------------------------- - * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] - * - * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity - * vector, and n the flow law exponent. - * - * If epsilon is NULL, it means this is the first time Emg is being run, and we - * return g, initial viscosity. - */ - - /*output: */ - IssmDouble viscosity3d; - - /*input strain rate: */ - IssmDouble exx,eyy,exy,exz,eyz; - - /*Intermediaries: */ - IssmDouble A,e; - IssmDouble B,n,D; - - /*Get B, Z and n*/ - n=GetN(); - D=GetD(); - B=(1-D)*GetB(); - - if (n==1){ - /*Viscous behaviour! viscosity3d=B: */ - viscosity3d=B/2; - } - else{ - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && - (epsilon[3]==0) && (epsilon[4]==0)){ - viscosity3d=0.5*pow(10.,14); - } - else{ - - /*Retrive strain rate components: */ - exx=*(epsilon+0); - eyy=*(epsilon+1); - exy=*(epsilon+2); - exz=*(epsilon+3); - eyz=*(epsilon+4); - - /*Build viscosity: viscosity3d=2*B/(2*A^e) */ - A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy; - if(A==0){ - /*Maxiviscosity3dm viscosity for 0 shear areas: */ - viscosity3d=2.25*pow(10.,17); - } - else{ - e=(n-1)/2/n; - - viscosity3d=B/(2*pow(A,e)); - } - } - } - - /*Checks in debugging mode*/ - if(viscosity3d<=0) _error_("Negative viscosity"); - _assert_(B>0); - _assert_(n>0); - - /*Assign output pointers:*/ - *pviscosity3d=viscosity3d; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosity3dFS {{{*/ -void Matdamageice::GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){ - /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: - * - * (1-D)* B - * viscosity3d= ------------------------------------------------------------------- - * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] - * - * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity - * vector, and n the flow law exponent. - * - * If epsilon is NULL, it means this is the first time Emg is being run, and we - * return g, initial viscosity. - */ - - /*output: */ - IssmDouble viscosity3d; - - /*input strain rate: */ - IssmDouble exx,eyy,exy,exz,eyz,ezz; - - /*Intermediaries: */ - IssmDouble A,e; - IssmDouble B,n,D; - IssmDouble eps0; - - /*Get B and n*/ - eps0=pow(10.,(IssmDouble)-27); - n=GetN(); - D=GetD(); - B=(1-D)*GetB(); - - if (n==1){ - /*Viscous behaviour! viscosity3d=B: */ - viscosity3d=B/2; - } - else{ - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && - (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ - viscosity3d=0.5*pow(10.,14); - } - else{ - - /*Retrive strain rate components: */ - exx=*(epsilon+0); - eyy=*(epsilon+1); - ezz=*(epsilon+2); //not used - exy=*(epsilon+3); - exz=*(epsilon+4); - eyz=*(epsilon+5); - - /*Build viscosity: viscosity3d=B/(2*A^e) */ - A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2); - if(A==0){ - /*Maxiviscosity3dm viscosity for 0 shear areas: */ - viscosity3d=2.25*pow(10.,17); - } - else{ - e=(n-1)/2/n; - viscosity3d=B/(2*pow(A,e)); - } - } - } - - /*Checks in debugging mode*/ - if(viscosity3d<=0) _error_("Negative viscosity"); - _assert_(B>0); - _assert_(n>0); - - /*Assign output pointers:*/ - *pviscosity3d=viscosity3d; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosityComplement {{{*/ -void Matdamageice::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ - /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: - * - * 1 - * viscosity= ------------------------------------------------------------------- - * 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(10.,17); - } - else{ - e=(n-1)/(2*n); - - viscosity_complement=1/(2*pow(A,e)); - } - } - else{ - viscosity_complement=4.5*pow(10.,17); - } - - /*Checks in debugging mode*/ - _assert_(B>0); - _assert_(n>0); - _assert_(viscosity_complement>0); - - /*Return: */ - *pviscosity_complement=viscosity_complement; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosityDComplement {{{*/ -void Matdamageice::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ - /*Return viscosity derivative for control method d(mu)/dD: - * - * 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(10.,17); - } - else{ - e=(n-1)/(2*n); - - viscosity_complement=- B/(2*pow(A,e)); - } - } - else{ - viscosity_complement=- 4.5*pow(10.,17); - } - - /*Checks in debugging mode*/ - _assert_(B>0); - _assert_(n>0); - _assert_(viscosity_complement<0); - - /*Return: */ - *pviscosity_complement=viscosity_complement; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosityDerivativeEpsSquare{{{*/ -void Matdamageice::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){ - - /*output: */ - IssmDouble mu_prime; - IssmDouble mu,n,eff2; - - /*input strain rate: */ - IssmDouble exx,eyy,exy,exz,eyz; - - /*Get visocisty and n*/ - GetViscosity3d(&mu,epsilon); - n=GetN(); - - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && - (epsilon[3]==0) && (epsilon[4]==0)){ - mu_prime=0.5*pow(10.,14); - } - else{ - /*Retrive strain rate components: */ - exx=epsilon[0]; - eyy=epsilon[1]; - exy=epsilon[2]; - exz=epsilon[3]; - eyz=epsilon[4]; - eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; - - mu_prime=(1-n)/(2*n) * mu/eff2; - } - - /*Assign output pointers:*/ - *pmu_prime=mu_prime; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosityDerivativeEpsSquareFS{{{*/ -void Matdamageice::GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* epsilon){ - - /*output: */ - IssmDouble mu_prime; - IssmDouble mu,n,eff2; - - /*input strain rate: */ - IssmDouble exx,eyy,ezz,exy,exz,eyz; - - /*Get visocisty and n*/ - GetViscosity3d(&mu,epsilon); - n=GetN(); - - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && - (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ - mu_prime=0.5*pow(10.,14); - } - else{ - /*Retrive strain rate components: */ - exx=epsilon[0]; - eyy=epsilon[1]; - ezz=epsilon[2]; - exy=epsilon[3]; - exz=epsilon[4]; - eyz=epsilon[5]; - eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; - - mu_prime=(1-n)/(2*n) * mu/eff2; - } - - /*Assign output pointers:*/ - *pmu_prime=mu_prime; -} -/*}}}*/ -/*FUNCTION Matdamageice::GetViscosity2dDerivativeEpsSquare{{{*/ -void Matdamageice::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){ - - /*output: */ - IssmDouble mu_prime; - IssmDouble mu,n,eff2; - - /*input strain rate: */ - IssmDouble exx,eyy,exy; - - /*Get visocisty and n*/ - GetViscosity2d(&mu,epsilon); - n=GetN(); - - if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ - mu_prime=0.5*pow(10.,14); - } - else{ - /*Retrive strain rate components: */ - exx=epsilon[0]; - eyy=epsilon[1]; - exy=epsilon[2]; - eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ; - - mu_prime=(1-n)/(2*n) * mu/eff2; - } - - /*Assign output pointers:*/ - *pmu_prime=mu_prime; -} -/*}}}*/ -/*FUNCTION Matdamageice::InputDuplicate{{{*/ -void Matdamageice::InputDuplicate(int original_enum,int new_enum){ - - /*Call inputs method*/ - if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); - -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/ -void Matdamageice::InputUpdateFromVector(IssmDouble* vector, int name, int type){ - - /*Intermediaries*/ - Element *element = NULL; - - /*Recover element*/ - element=(Element*)helement->delivers(); - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - switch(type){ - - case VertexEnum: - - switch(element->ObjectEnum()){ - - case TriaEnum: { - IssmDouble values[3]; - for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Pid()]; - this->inputs->AddInput(new TriaInput(name,values,P1Enum)); - return; - } - default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet"); - } - default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); - } -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/ -void Matdamageice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){ - - /*Intermediaries*/ - Element *element = NULL; - Parameters* parameters= NULL; - int dim; - - /*Recover element*/ - element=(Element*)helement->delivers(); - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - switch(type){ - - case VertexEnum: - - switch(element->ObjectEnum()){ - - case TriaEnum: { - IssmDouble values[3]; - for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector - this->inputs->AddInput(new TriaInput(name,values,P1Enum)); - /*Special case for rheology B in 2D: Pourave land for this solution{{{*/ - if(name==MaterialsRheologyBEnum){ - /*Are we in 2D?:*/ - if(element->ObjectEnum()==TriaEnum){ - parameters=((Tria*)(element))->parameters; - } - else{ - parameters=((Penta*)(element))->parameters; - } - parameters->FindParam(&dim,MeshDimensionEnum); - if(dim==2){ - /*Dupliacte rheology input: */ - this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values,P1Enum)); - } - } - /*}}}*/ - return; - } - default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet"); - } - default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); - } - -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromMatrixDakota(IssmDouble* vector, int name, int type) {{{*/ -void Matdamageice::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/ -void Matdamageice::InputUpdateFromConstant(IssmDouble constant, int name){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromConstant(int constant, int name) {{{*/ -void Matdamageice::InputUpdateFromConstant(int constant, int name){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromConstant(bool constant, int name) {{{*/ -void Matdamageice::InputUpdateFromConstant(bool constant, int name){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromSolution{{{*/ -void Matdamageice::InputUpdateFromSolution(IssmDouble* solution){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Matdamageice::InputUpdateFromIoModel{{{*/ -void Matdamageice::InputUpdateFromIoModel(int index, IoModel* iomodel){ - - int i,j; - bool control_analysis; - int num_control_type; - - /*Fetch parameters: */ - iomodel->Constant(&control_analysis,InversionIscontrolEnum); - if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum); - - if(iomodel->dim==2){ - - /*Intermediaries*/ - const int num_vertices = 3; //Tria has 3 vertices - IssmDouble nodeinputs[num_vertices]; - IssmDouble cmmininputs[num_vertices]; - IssmDouble cmmaxinputs[num_vertices]; - - /*Get B*/ - if (iomodel->Data(MaterialsRheologyBEnum)) { - for(i=0;iData(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+i]-1]; - this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,nodeinputs,P1Enum)); - } - - /*Get n*/ - if (iomodel->Data(MaterialsRheologyNEnum)) { - for(i=0;iData(MaterialsRheologyNEnum)[index]; - this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum)); - } - - /*Get D*/ - if (iomodel->Data(DamageDEnum)) { - for(i=0;iData(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1]; - this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum)); - } - - /*Control Inputs*/ - #ifdef _HAVE_CONTROL_ - if (control_analysis && iomodel->Data(InversionControlParametersEnum)){ - for(i=0;i(iomodel->Data(InversionControlParametersEnum)[i])){ - case MaterialsRheologyBbarEnum: - if (iomodel->Data(MaterialsRheologyBEnum)){ - _assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); - for(j=0;jData(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+j]-1]; - for(j=0;jData(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - for(j=0;jData(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); - } - break; - case DamageDbarEnum: - if (iomodel->Data(DamageDEnum)){ - _assert_(iomodel->Data(DamageDEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); - for(j=0;jData(DamageDEnum)[iomodel->elements[num_vertices*index+j]-1]; - for(j=0;jData(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - for(j=0;jData(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); - } - break; - - } - } - } - #endif - } - - /*if 3d*/ - #ifdef _HAVE_3D_ - else if(iomodel->dim==3){ - - /*Intermediaries*/ - const int num_vertices = 6; //Penta has 6 vertices - IssmDouble nodeinputs[num_vertices]; - IssmDouble cmmininputs[num_vertices]; - IssmDouble cmmaxinputs[num_vertices]; - - /*Get B*/ - if (iomodel->Data(MaterialsRheologyBEnum)) { - for(i=0;iData(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+i]-1]; - this->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum)); - } - - /*Get n*/ - if (iomodel->Data(MaterialsRheologyNEnum)) { - for(i=0;iData(MaterialsRheologyNEnum)[index]; - this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum)); - } - - /*Get D*/ - if (iomodel->Data(DamageDEnum)) { - for(i=0;iData(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1]; - this->inputs->AddInput(new PentaInput(DamageDEnum,nodeinputs,P1Enum)); - } - - /*Control Inputs*/ - #ifdef _HAVE_CONTROL_ - if (control_analysis && iomodel->Data(InversionControlParametersEnum)){ - for(i=0;i(iomodel->Data(InversionControlParametersEnum)[i])){ - case MaterialsRheologyBbarEnum: - if (iomodel->Data(MaterialsRheologyBEnum)){ - _assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); - for(j=0;jData(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+j]-1]; - for(j=0;jData(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - for(j=0;jData(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); - } - break; - case DamageDbarEnum: - if (iomodel->Data(DamageDEnum)){ - _assert_(iomodel->Data(DamageDEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); - for(j=0;jData(DamageDEnum)[iomodel->elements[num_vertices*index+j]-1]; - for(j=0;jData(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - for(j=0;jData(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); - } - break; - } - } - } - #endif - } - #endif - else{ - _error_("Mesh type not supported yet!"); - } - - return; -} -/*}}}*/ -/*FUNCTION Matdamageice::IsInput{{{*/ -bool Matdamageice::IsInput(int name){ - if ( - name==MaterialsRheologyBEnum || - name==MaterialsRheologyBbarEnum || - name==MaterialsRheologyNEnum || - name==DamageDEnum || - name==DamageDbarEnum - ){ - return true; - } - else return false; -} -/*}}}*/ Index: sm/trunk-jpl/src/c/classes/Materials/Matdamageice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.h (revision 16166) +++ (revision ) @@ -1,68 +1,0 @@ -/*!\file Matdamageice.h - * \brief: header file for matice object - */ - -#ifndef MATDAMAGEICE_H_ -#define MATDAMAGEICE_H_ - -/*Headers:*/ -/*{{{*/ -#include "./Material.h" -class IoModel; -/*}}}*/ - -class Matdamageice: public Material{ - - private: - int mid; - Hook* helement; - - public: - /*Matdamageice constructors, destructors: {{{*/ - Matdamageice(); - Matdamageice(int mid,int i, IoModel* iomodel); - ~Matdamageice(); - /*}}}*/ - /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); - int ObjectEnum(); - Object* copy(); - /*}}}*/ - /*Update virtual functions definitions: {{{*/ - void InputUpdateFromVector(IssmDouble* vector, int name, int type); - void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols, int name, int type); - void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); - void InputUpdateFromConstant(IssmDouble constant, int name); - void InputUpdateFromConstant(int constant, int name); - void InputUpdateFromConstant(bool constant, int name); - void InputUpdateFromSolution(IssmDouble* solution); - void InputUpdateFromIoModel(int index, IoModel* iomodel); - /*}}}*/ - /*Material virtual functions resolution: {{{*/ - void InputDuplicate(int original_enum,int new_enum); - void Configure(Elements* elements); - void GetVectorFromInputs(Vector* vector,int input_enum); - 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 GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon); - void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); - void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon); - void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); - IssmDouble GetA(); - IssmDouble GetB(); - IssmDouble GetBbar(); - IssmDouble GetN(); - IssmDouble GetZ(); - IssmDouble GetD(); - IssmDouble GetZbar(); - IssmDouble GetDbar(); - bool IsInput(int name); - /*}}}*/ -}; - -#endif /* _MATICE_H_ */ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16167) @@ -175,4 +175,42 @@ } /*}}}*/ +/*FUNCTION Matice::GetD {{{*/ +IssmDouble Matice::GetD(){ + + /*Output*/ + IssmDouble D; + + inputs->GetInputAverage(&D,DamageDEnum); + return D; +} +/*}}}*/ +/*FUNCTION Matice::GetZ {{{*/ +IssmDouble Matice::GetZ(){ + + /*Output*/ + IssmDouble D; + + inputs->GetInputAverage(&D,DamageDEnum); + return 1/(1-D); +} +/*}}}*/ +/*FUNCTION Matice::GetZbar {{{*/ +IssmDouble Matice::GetZbar(){ + + /*Output*/ + IssmDouble Dbar; + inputs->GetInputAverage(&Dbar,DamageDbarEnum); + return 1/(1-Dbar); +} +/*}}}*/ +/*FUNCTION Matice::GetDbar {{{*/ +IssmDouble Matice::GetDbar(){ + + /*Output*/ + IssmDouble Dbar; + inputs->GetInputAverage(&Dbar,DamageDbarEnum); + return Dbar; +} +/*}}}*/ /*FUNCTION Matice::GetVectorFromInputs{{{*/ void Matice::GetVectorFromInputs(Vector* vector,int input_enum){ @@ -210,5 +248,5 @@ void Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ /*From a string tensor and a material object, return viscosity, using Glen's flow law. - B + (1-D) B viscosity= ------------------------------------------------------------------- 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -229,13 +267,14 @@ /*Intermediary: */ IssmDouble A,e; - IssmDouble B,n; + IssmDouble B,D,n; /*Get B and n*/ B=GetBbar(); n=GetN(); + D=GetDbar(); if (n==1){ /*Viscous behaviour! viscosity=B: */ - viscosity=B/2; + viscosity=(1-D)*B/2; } else{ @@ -257,5 +296,5 @@ else{ e=(n-1)/(2*n); - viscosity=B/(2*pow(A,e)); + viscosity=(1-D)*B/(2*pow(A,e)); } } @@ -266,4 +305,5 @@ _assert_(B>0); _assert_(n>0); + _assert_(D>=0 && D<1); /*Return: */ @@ -276,5 +316,5 @@ /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: * - * B + * (1-D)*B * viscosity3d= ------------------------------------------------------------------- * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -295,13 +335,14 @@ /*Intermediaries: */ IssmDouble A,e; - IssmDouble B,n; + IssmDouble B,D,n; /*Get B and n*/ B=GetB(); + D=GetD(); n=GetN(); if (n==1){ /*Viscous behaviour! viscosity3d=B: */ - viscosity3d=B/2; + viscosity3d=(1-D)*B/2; } else{ @@ -328,5 +369,5 @@ e=(n-1)/2/n; - viscosity3d=B/(2*pow(A,e)); + viscosity3d=(1-D)*B/(2*pow(A,e)); } } @@ -337,4 +378,5 @@ _assert_(B>0); _assert_(n>0); + _assert_(D>=0 && D<1); /*Assign output pointers:*/ @@ -346,5 +388,5 @@ /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: * - * B + * (1-D)*B * viscosity3d= ------------------------------------------------------------------- * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -365,5 +407,5 @@ /*Intermediaries: */ IssmDouble A,e; - IssmDouble B,n; + IssmDouble B,D,n; IssmDouble eps0; @@ -371,9 +413,10 @@ eps0=pow(10.,-27); B=GetB(); + D=GetD(); n=GetN(); if (n==1){ /*Viscous behaviour! viscosity3d=B: */ - viscosity3d=B/2; + viscosity3d=(1-D)*B/2; } else{ @@ -400,5 +443,5 @@ else{ e=(n-1)/2/n; - viscosity3d=B/(2*pow(A,e)); + viscosity3d=(1-D)*B/(2*pow(A,e)); } } @@ -409,4 +452,5 @@ _assert_(B>0); _assert_(n>0); + _assert_(D>=0 && D<1); /*Assign output pointers:*/ @@ -418,5 +462,5 @@ /*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: * - * 1 + * (1-D) * viscosity= ------------------------------------------------------------------- * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -434,8 +478,8 @@ /*Intermediary value A and exponent e: */ IssmDouble A,e; - IssmDouble B,n; - - /*Get B and n*/ - B=GetBbar(); + IssmDouble D,n; + + /*Get D and n*/ + D=GetDbar(); n=GetN(); @@ -445,18 +489,74 @@ exy=*(epsilon+2); + /*Build viscosity: mu2=(1-D)/(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(10.,17); + } + else{ + e=(n-1)/(2*n); + + viscosity_complement=(1-D)/(2*pow(A,e)); + } + } + else{ + viscosity_complement=4.5*pow(10.,17); + } + + /*Checks in debugging mode*/ + _assert_(D>=0 && D<1); + _assert_(n>0); + _assert_(viscosity_complement>0); + + /*Return: */ + *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosityDComplement {{{*/ +void Matice::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ + /*Return viscosity derivative for control method d(mu)/dD: + * + * 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(10.,17); + viscosity_complement=- 2.25*pow(10.,17); } else{ e=(n-1)/(2*n); - viscosity_complement=1/(2*pow(A,e)); + viscosity_complement=- B/(2*pow(A,e)); } } else{ - viscosity_complement=4.5*pow(10.,17); + viscosity_complement=- 4.5*pow(10.,17); } @@ -464,5 +564,5 @@ _assert_(B>0); _assert_(n>0); - _assert_(viscosity_complement>0); + _assert_(viscosity_complement<0); /*Return: */ @@ -729,4 +829,10 @@ } + /*Get D:*/ + if (iomodel->Data(DamageDEnum)) { + for(i=0;iData(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1]; + this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum)); + } + /*Control Inputs*/ #ifdef _HAVE_CONTROL_ @@ -743,4 +849,13 @@ } break; + case DamageDbarEnum: + if (iomodel->Data(DamageDEnum)){ + _assert_(iomodel->Data(DamageDEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(DamageDEnum)[iomodel->elements[num_vertices*index+j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; } } @@ -770,5 +885,11 @@ this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum)); } - + + /*Get D*/ + if (iomodel->Data(DamageDEnum)) { + for(i=0;iData(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1]; + this->inputs->AddInput(new PentaInput(DamageDEnum,nodeinputs,P1Enum)); + } + /*Control Inputs*/ #ifdef _HAVE_CONTROL_ @@ -785,4 +906,14 @@ } break; + case DamageDbarEnum: + if (iomodel->Data(DamageDEnum)){ + _assert_(iomodel->Data(DamageDEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(DamageDEnum)[iomodel->elements[num_vertices*index+j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + } } @@ -804,5 +935,8 @@ name==MaterialsRheologyBbarEnum || name==MaterialsRheologyNEnum || - name==QmuMaterialsRheologyBEnum + name==QmuMaterialsRheologyBEnum || + name==DamageDEnum || + name==DamageDbarEnum + ){ return true; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 16166) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 16167) @@ -57,5 +57,5 @@ void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon); void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityDComplement(IssmDouble*, IssmDouble*){_error_("not supported");}; + void GetViscosityDComplement(IssmDouble*, IssmDouble*); void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon); @@ -64,8 +64,8 @@ IssmDouble GetB(); IssmDouble GetBbar(); - IssmDouble GetZ(){_error_("not supported");}; - IssmDouble GetD(){_error_("not supported");}; - IssmDouble GetZbar(){_error_("not supported");}; - IssmDouble GetDbar(){_error_("not supported");}; + IssmDouble GetZ(); + IssmDouble GetD(); + IssmDouble GetZbar(); + IssmDouble GetDbar(); IssmDouble GetN(); bool IsInput(int name); Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 16166) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 16167) @@ -77,5 +77,4 @@ #include "./Materials/Material.h" #include "./Materials/Matice.h" -#include "./Materials/Matdamageice.h" #include "./Materials/Matpar.h" Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 16167) @@ -34,4 +34,5 @@ for(j=0;jSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); element->Gradj(gradient_list[i],control_type[i],i); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 16167) @@ -54,5 +54,5 @@ switch(materials_type){ case MaticeEnum: - iomodel->FetchData(2,MaterialsRheologyBEnum,MaterialsRheologyNEnum); + iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,DamageDEnum); for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); if(dakota_analysis){ @@ -62,8 +62,4 @@ #endif } - break; - case MatdamageiceEnum: - iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,DamageDEnum); - for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matdamageice(i+1,i,iomodel)); break; default: Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp (revision 16167) @@ -101,7 +101,6 @@ iomodel->FetchDataToInput(elements,LoadingforceXEnum); iomodel->FetchDataToInput(elements,LoadingforceYEnum); - if(materials_type==MatdamageiceEnum){ - iomodel->FetchDataToInput(elements,DamageDEnum); - } + iomodel->FetchDataToInput(elements,DamageDEnum); + if(iomodel->dim==3){ iomodel->FetchDataToInput(elements,MeshElementonbedEnum); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16166) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16167) @@ -152,6 +152,4 @@ MaterialsRheologyLawEnum, MaterialsRheologyNEnum, - MaterialsRheologyZEnum, - MaterialsRheologyZbarEnum, DamageDEnum, DamageDbarEnum, @@ -373,5 +371,4 @@ SSA3dIceFrontEnum, MaticeEnum, - MatdamageiceEnum, MatparEnum, NodeEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16167) @@ -160,6 +160,4 @@ case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; case MaterialsRheologyNEnum : return "MaterialsRheologyN"; - case MaterialsRheologyZEnum : return "MaterialsRheologyZ"; - case MaterialsRheologyZbarEnum : return "MaterialsRheologyZbar"; case DamageDEnum : return "DamageD"; case DamageDbarEnum : return "DamageDbar"; @@ -373,5 +371,4 @@ case SSA3dIceFrontEnum : return "SSA3dIceFront"; case MaticeEnum : return "Matice"; - case MatdamageiceEnum : return "Matdamageice"; case MatparEnum : return "Matpar"; case NodeEnum : return "Node"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16166) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16167) @@ -163,6 +163,4 @@ 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,"DamageD")==0) return DamageDEnum; else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum; @@ -260,10 +258,10 @@ else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum; else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; + else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum; + else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum; - else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum; - else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; + if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum; else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; @@ -379,15 +377,14 @@ else if (strcmp(name,"SSA3dIceFront")==0) return SSA3dIceFrontEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; - else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matpar")==0) return MatparEnum; else if (strcmp(name,"Node")==0) return NodeEnum; else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; + else if (strcmp(name,"NumericalfluxType")==0) return NumericalfluxTypeEnum; + else if (strcmp(name,"Param")==0) return ParamEnum; + else if (strcmp(name,"L1L2IceFront")==0) return L1L2IceFrontEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"NumericalfluxType")==0) return NumericalfluxTypeEnum; - else if (strcmp(name,"Param")==0) return ParamEnum; - else if (strcmp(name,"L1L2IceFront")==0) return L1L2IceFrontEnum; - else if (strcmp(name,"HOIceFront")==0) return HOIceFrontEnum; + if (strcmp(name,"HOIceFront")==0) return HOIceFrontEnum; else if (strcmp(name,"Pengrid")==0) return PengridEnum; else if (strcmp(name,"Penpair")==0) return PenpairEnum; @@ -506,11 +503,11 @@ else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; + else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; + else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; + else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; - else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; - else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; - else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; + if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; else if (strcmp(name,"GiaW")==0) return GiaWEnum; else if (strcmp(name,"P0")==0) return P0Enum; Index: /issm/trunk-jpl/src/m/classes/damage.m =================================================================== --- /issm/trunk-jpl/src/m/classes/damage.m (revision 16166) +++ /issm/trunk-jpl/src/m/classes/damage.m (revision 16167) @@ -56,6 +56,6 @@ function md = checkconsistency(obj,md,solution,analyses) % {{{ + md = checkfield(md,'damage.D','>=0',0,'size',[md.mesh.numberofvertices 1]); if obj.isdamage, - md = checkfield(md,'damage.D','>=0',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'damage.law','values',{'undamaged','pralong'}); end @@ -86,6 +86,6 @@ WriteData(fid,'object',obj,'class','damage','fieldname','isdamage','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','damage','fieldname','D','format','DoubleMat','mattype',1); if obj.isdamage, - WriteData(fid,'object',obj,'class','damage','fieldname','D','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','damage','fieldname','law','format','String'); end Index: /issm/trunk-jpl/src/m/classes/inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.m (revision 16166) +++ /issm/trunk-jpl/src/m/classes/inversion.m (revision 16167) @@ -116,5 +116,5 @@ md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]); md = checkfield(md,'inversion.control_parameters','cell',1,'values',... - {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsDamageDbar' 'Vx' 'Vy' 'Thickness'}); + {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness'}); md = checkfield(md,'inversion.nsteps','numel',1,'>=',1); md = checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0); Index: /issm/trunk-jpl/src/m/classes/inversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.py (revision 16166) +++ /issm/trunk-jpl/src/m/classes/inversion.py (revision 16167) @@ -124,5 +124,5 @@ 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','MaterialsDamageDbar','Vx','Vy']) + md = checkfield(md,'inversion.control_parameters','cell',1,'values',['BalancethicknessThickeningRate','FrictionCoefficient','MaterialsRheologyBbar','DamageDbar','Vx','Vy']) md = checkfield(md,'inversion.nsteps','numel',[1],'>=',1) md = checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps],'>=',0) Index: sm/trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 16166) +++ (revision ) @@ -1,209 +1,0 @@ -%MATDAMAGEICE class definition -% -% Usage: -% matdamageice=matdamageice(); - -classdef matdamageice - properties (SetAccess=public) - % {{{ - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_B = NaN; - rheology_n = NaN; - rheology_law = ''; - - %damage - damage_D = NaN; - damage_law = ''; - - %parameters for law 'initial': - damage_stress_threshold = NaN; - damage_c1 = NaN; - damage_c2 = NaN; - damage_c3 = NaN; - damage_c4 = NaN; - - %gia: - lithosphere_shear_modulus = 0.; - lithosphere_density = 0.; - mantle_shear_modulus = 0.; - mantle_density = 0.; - - end % }}} - methods - function obj = matdamageice(varargin) % {{{ - switch nargin - case 0 - obj=setdefaultparameters(obj); - case 1 - inputstruct=varargin{1}; - list1 = properties('matdamageice'); - list2 = fieldnames(inputstruct); - for i=1:length(list1) - fieldname = list1{i}; - if ismember(fieldname,list2), - obj.(fieldname) = inputstruct.(fieldname); - end - end - otherwise - error('constructor not supported'); - end - end % }}} - function obj = setdefaultparameters(obj) % {{{ - - %ice density (kg/m^3) - obj.rho_ice=917.; - - %ocean water density (kg/m^3) - obj.rho_water=1023.; - - %fresh water density (kg/m^3) - obj.rho_freshwater=1000.; - - %water viscosity (N.s/m^2) - obj.mu_water=0.001787; - - %ice heat capacity cp (J/kg/K) - obj.heatcapacity=2093.; - - %ice latent heat of fusion L (J/kg) - obj.latentheat=3.34*10^5; - - %ice thermal conductivity (W/m/K) - obj.thermalconductivity=2.4; - - %the melting point of ice at 1 atmosphere of pressure in K - obj.meltingpoint=273.15; - - %rate of change of melting point with pressure (K/Pa) - obj.beta=9.8*10^-8; - - %mixed layer (ice-water interface) heat capacity (J/kg/K) - obj.mixed_layer_capacity=3974.; - - %thermal exchange velocity (ice-water interface) (m/s) - obj.thermal_exchange_velocity=1.00*10^-4; - - %Rheology law: what is the temperature dependence of B with T - %available: none, paterson and arrhenius - obj.rheology_law='Paterson'; - - %damage parameters: - obj.damage_law='initial'; - obj.damage_stress_threshold=0; - obj.damage_c1=0; - obj.damage_c2=0; - obj.damage_c3=0; - obj.damage_c4=0; - - %GIA: - obj.lithosphere_shear_modulus = 6.7*10^10; % (Pa) - obj.lithosphere_density = 3.32; % (g/cm^-3) - obj.mantle_shear_modulus = 1.45*10^11; % (Pa) - obj.mantle_density = 3.34; % (g/cm^-3) - - end % }}} - function md = checkconsistency(obj,md,solution,analyses) % {{{ - md = checkfield(md,'materials.rho_ice','>',0); - md = checkfield(md,'materials.rho_water','>',0); - md = checkfield(md,'materials.rho_freshwater','>',0); - md = checkfield(md,'materials.mu_water','>',0); - 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_law','values',{'None' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); - - md = checkfield(md,'damage.D','>',0,'size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'damage.law','values',{'undamaged' 'pralong'}); - if strcmpi(obj.damage_law,'initial'), - md = checkfield(md,'damage.c1','>=',0); - md = checkfield(md,'damage.c2','>=',0); - md = checkfield(md,'damage.c3','>=',0); - md = checkfield(md,'damage.c4','>=',0); - md = checkfield(md,'damage.stress_threshold','>=',0); - end - - if ismember(GiaAnalysisEnum(),analyses), - md = checkfield(md,'materials.lithosphere_shear_modulus','>',0,'numel',1); - md = checkfield(md,'materials.lithosphere_density','>',0,'numel',1); - md = checkfield(md,'materials.mantle_shear_modulus','>',0,'numel',1); - md = checkfield(md,'materials.mantle_density','>',0,'numel',1); - end - - end % }}} - function disp(obj) % {{{ - disp(sprintf(' Materials:\n')); - - fielddisplay(obj,'rho_ice','ice density [kg/m^3]'); - fielddisplay(obj,'rho_water','ocean water density [kg/m^3]'); - fielddisplay(obj,'rho_freshwater','fresh water density [kg/m^3]'); - fielddisplay(obj,'mu_water','water viscosity [N s/m^2]'); - fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]'); - fielddisplay(obj,'thermalconductivity','ice thermal conductivity [W/m/K]'); - fielddisplay(obj,'meltingpoint','melting point of ice at 1atm in K'); - fielddisplay(obj,'latentheat','latent heat of fusion [J/m^3]'); - fielddisplay(obj,'beta','rate of change of melting point with pressure [K/Pa]'); - fielddisplay(obj,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); - fielddisplay(obj,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); - fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]'); - fielddisplay(obj,'rheology_n','Glen''s flow law exponent'); - fielddisplay(obj,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); - fielddisplay(obj,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); - fielddisplay(obj,'lithosphere_density','Lithosphere density [g/cm^-3]'); - fielddisplay(obj,'mantle_shear_modulus','Mantle shear modulus [Pa]'); - fielddisplay(obj,'mantle_density','Mantle density [g/cm^-3]'); - fielddisplay(obj,'D','damage tensor (scalar)'); - fielddisplay(obj,'law','damage law (string) from {''initial''}'); - if strcmpi(obj.damage_law,'initial'), - fielddisplay(obj,'c1','damage parameter 1'); - fielddisplay(obj,'c2','damage parameter 2'); - fielddisplay(obj,'c3','damage parameter 3'); - fielddisplay(obj,'c4','damage parameter 4'); - fielddisplay(obj,'stress_threshold','damage stress threshold [Pa]'); - end - - end % }}} - function marshall(obj,md,fid) % {{{ - WriteData(fid,'enum',MaterialsEnum(),'data',MatdamageiceEnum(),'format','Integer'); - WriteData(fid,'object',obj,'class','materials','fieldname','rho_ice','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','rho_water','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','rho_freshwater','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','mu_water','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','heatcapacity','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','latentheat','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','thermalconductivity','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','meltingpoint','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','beta','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','mixed_layer_capacity','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); - WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer'); - - WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_shear_modulus','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3); - WriteData(fid,'object',obj,'class','materials','fieldname','mantle_shear_modulus','format','Double'); - WriteData(fid,'object',obj,'class','materials','fieldname','mantle_density','format','Double','scale',10^3); - - %damage: - WriteData(fid,'object',obj,'class','damage','fieldname','D','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','damage','fieldname','law','format','String'); - if strcmpi(obj.damage_law,'pralong'), - WriteData(fid,'object',obj,'class','damage','fieldname','c1','format','Double'); - WriteData(fid,'object',obj,'class','damage','fieldname','c2','format','Double'); - WriteData(fid,'object',obj,'class','damage','fieldname','c3','format','Double'); - WriteData(fid,'object',obj,'class','damage','fieldname','c4','format','Double'); - WriteData(fid,'object',obj,'class','damage','fieldname','stress_threshold','format','Double'); - end - - end % }}} - end -end Index: sm/trunk-jpl/src/m/classes/matdamageice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 16166) +++ (revision ) @@ -1,193 +1,0 @@ -from fielddisplay import fielddisplay -from EnumDefinitions import * -from StringToEnum import StringToEnum -from checkfield import * -from WriteData import * - -class matdamageice(object): - """ - MATDAMAGEICE class definition - - Usage: - matdamageice=matdamageice(); - """ - - def __init__(self,*args): # {{{ - self.rho_ice = 0. - self.rho_water = 0. - self.rho_freshwater = 0. - self.mu_water = 0. - self.heatcapacity = 0. - self.latentheat = 0. - self.thermalconductivity = 0. - self.meltingpoint = 0. - self.beta = 0. - self.mixed_layer_capacity = 0. - self.thermal_exchange_velocity = 0. - self.rheology_B = float('NaN') - self.rheology_n = float('NaN') - self.rheology_law = '' - - #damage: - self.damage_D = float('NaN') - self.damage_law = '' - - #parameters for law 'initial': - self.damage_c1 = float('NaN') - self.damage_c2 = float('NaN') - self.damage_c3 = float('NaN') - self.damage_c4 = float('NaN') - self.damage_stress_threshold = float('NaN') - - #gia: - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. - - if not len(args): - self.setdefaultparameters() - else: - raise RuntimeError("constructor not supported") - - # }}} - def __repr__(self): # {{{ - s =' Materials:\n' - - s+="%s\n" % fielddisplay(self,"rho_ice","ice density [kg/m^3]") - s+="%s\n" % fielddisplay(self,"rho_water","ocean water density [kg/m^3]") - s+="%s\n" % fielddisplay(self,"rho_freshwater","fresh water density [kg/m^3]") - s+="%s\n" % fielddisplay(self,"mu_water","water viscosity [N s/m^2]") - s+="%s\n" % fielddisplay(self,"heatcapacity","heat capacity [J/kg/K]") - s+="%s\n" % fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]") - s+="%s\n" % fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K") - s+="%s\n" % fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]") - s+="%s\n" % fielddisplay(self,"beta","rate of change of melting point with pressure [K/Pa]") - s+="%s\n" % fielddisplay(self,"mixed_layer_capacity","mixed layer capacity [W/kg/K]") - s+="%s\n" % fielddisplay(self,"thermal_exchange_velocity","thermal exchange velocity [m/s]") - s+="%s\n" % fielddisplay(self,"rheology_B","flow law parameter [Pa/s^(1/n)]") - s+="%s\n" % fielddisplay(self,"rheology_n","Glen's flow law exponent") - s+="%s\n" % fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'Paterson', 'Arrhenius' or 'LliboutryDuval'") - s+="%s\n" % fielddisplay(self,"lithosphere_shear_modulus","Lithosphere shear modulus [Pa]") - s+="%s\n" % fielddisplay(self,"lithosphere_density","Lithosphere density [g/cm^-3]") - s+="%s\n" % fielddisplay(self,"mantle_shear_modulus","Mantle shear modulus [Pa]") - s+="%s\n" % fielddisplay(self,"mantle_density","Mantle density [g/cm^-3]") - s+="%s\n" % fielddisplay(self,"damage_D","damage tensor (scalar for now)") - - if (obj.damage_law=='initial'): - s+="%s\n" % fielddisplay(self,"damage_c1","damage parameter 1 ") - s+="%s\n" % fielddisplay(self,"damage_c2","damage parameter 2 ") - s+="%s\n" % fielddisplay(self,"damage_c3","damage parameter 3 ") - s+="%s\n" % fielddisplay(self,"damage_c4","damage parameter 4 ") - s+="%s\n" % fielddisplay(self,"damage_stress_threshold","damage stress threshold [Pa]") - - return s - # }}} - def setdefaultparameters(self): # {{{ - - #ice density (kg/m^3) - self.rho_ice=917. - - #ocean water density (kg/m^3) - self.rho_water=1023. - - #fresh water density (kg/m^3) - self.rho_freshwater=1000. - - #water viscosity (N.s/m^2) - self.mu_water=0.001787 - - #ice heat capacity cp (J/kg/K) - self.heatcapacity=2093. - - #ice latent heat of fusion L (J/kg) - self.latentheat=3.34*10**5 - - #ice thermal conductivity (W/m/K) - self.thermalconductivity=2.4 - - #the melting point of ice at 1 atmosphere of pressure in K - self.meltingpoint=273.15 - - #rate of change of melting point with pressure (K/Pa) - self.beta=9.8*10**-8 - - #mixed layer (ice-water interface) heat capacity (J/kg/K) - self.mixed_layer_capacity=3974. - - #thermal exchange velocity (ice-water interface) (m/s) - self.thermal_exchange_velocity=1.00*10**-4 - - #Rheology law: what is the temperature dependence of B with T - #available: none, paterson and arrhenius - self.rheology_law='Paterson' - - #damage parameters: - self.damage_law='initial' - self.damage_stress_threshold=0 - self.damage_c1=0 - self.damage_c2=0 - self.damage_c3=0 - self.damage_c4=0 - - # GIA: - self.lithosphere_shear_modulus = 6.7*10**10; # (Pa) - self.lithosphere_density = 3.32; # (g/cm^-3) - self.mantle_shear_modulus = 1.45*10**11; # (Pa) - self.mantle_density = 3.34; # (g/cm^-3) - # }}} - def checkconsistency(self,md,solution,analyses): # {{{ - md = checkfield(md,'materials.rho_ice','>',0) - md = checkfield(md,'materials.rho_water','>',0) - md = checkfield(md,'materials.rho_freshwater','>',0) - md = checkfield(md,'materials.mu_water','>',0) - md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements]) - md = checkfield(md,'materials.rheology_law','values',['None','Paterson','Arrhenius','LliboutryDuval']) - md = checkfield(md,'materials.lithosphere_shear_modulus','>',0,'numel',[1]); - md = checkfield(md,'materials.lithosphere_density','>',0,'numel',[1]); - md = checkfield(md,'materials.mantle_shear_modulus','>',0,'numel',[1]); - md = checkfield(md,'materials.mantle_density','>',0,'numel',[1]); - - md = checkfield(md,'materials.damage_D','>',0,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'materials.damage_law','values',['initial']) - if self.damage_law == 'initial': - md = checkfield(md,'materials.damage_c1','>=',0) - md = checkfield(md,'materials.damage_c2','>=',0) - md = checkfield(md,'materials.damage_c3','>=',0) - md = checkfield(md,'materials.damage_c4','>=',0) - md = checkfield(md,'materials.damage_stress_threshold','>=',0) - - return md - # }}} - def marshall(self,md,fid): # {{{ - WriteData(fid,'enum',MaterialsEnum(),'data',MatdamageiceEnum(),'format','Integer') - WriteData(fid,'object',self,'class','materials','fieldname','rho_ice','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','rho_water','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','rho_freshwater','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','mu_water','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','heatcapacity','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','latentheat','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','thermalconductivity','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','meltingpoint','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','beta','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double') - WriteData(fid,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1) - WriteData(fid,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2) - WriteData(fid,'data',StringToEnum(self.rheology_law)[0],'enum',MaterialsRheologyLawEnum(),'format','Integer') - - WriteData(fid,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double'); - WriteData(fid,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10.**3.); - WriteData(fid,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double'); - WriteData(fid,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10.**3.); - - WriteData(fid,'object',self,'class','damage','fieldname','D','format','DoubleMat','mattype',1) - WriteData(fid,'object',self,'class','damage','fieldname','law','format','String') - if self.damage_law=='initial': - WriteData(fid,'object',self,'class','damage','fieldname','c1','format','Double') - WriteData(fid,'object',self,'class','damage','fieldname','c2','format','Double') - WriteData(fid,'object',self,'class','damage','fieldname','c3','format','Double') - WriteData(fid,'object',self,'class','damage','fieldname','c4','format','Double') - WriteData(fid,'object',self,'class','damage','fieldname','stress_threshold','format','Double') - # }}} Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16166) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16167) @@ -152,6 +152,4 @@ def MaterialsRheologyLawEnum(): return StringToEnum("MaterialsRheologyLaw")[0] def MaterialsRheologyNEnum(): return StringToEnum("MaterialsRheologyN")[0] -def MaterialsRheologyZEnum(): return StringToEnum("MaterialsRheologyZ")[0] -def MaterialsRheologyZbarEnum(): return StringToEnum("MaterialsRheologyZbar")[0] def DamageDEnum(): return StringToEnum("DamageD")[0] def DamageDbarEnum(): return StringToEnum("DamageDbar")[0] @@ -365,5 +363,4 @@ def SSA3dIceFrontEnum(): return StringToEnum("SSA3dIceFront")[0] def MaticeEnum(): return StringToEnum("Matice")[0] -def MatdamageiceEnum(): return StringToEnum("Matdamageice")[0] def MatparEnum(): return StringToEnum("Matpar")[0] def NodeEnum(): return StringToEnum("Node")[0] Index: sm/trunk-jpl/src/m/enum/MatdamageiceEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MatdamageiceEnum.m (revision 16166) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=MatdamageiceEnum() -%MATDAMAGEICEENUM - Enum of Matdamageice -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=MatdamageiceEnum() - -macro=StringToEnum('Matdamageice'); Index: sm/trunk-jpl/src/m/enum/MaterialsRheologyZEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaterialsRheologyZEnum.m (revision 16166) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=MaterialsRheologyZEnum() -%MATERIALSRHEOLOGYZENUM - Enum of MaterialsRheologyZ -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=MaterialsRheologyZEnum() - -macro=StringToEnum('MaterialsRheologyZ'); Index: sm/trunk-jpl/src/m/enum/MaterialsRheologyZbarEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaterialsRheologyZbarEnum.m (revision 16166) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=MaterialsRheologyZbarEnum() -%MATERIALSRHEOLOGYZBARENUM - Enum of MaterialsRheologyZbar -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=MaterialsRheologyZbarEnum() - -macro=StringToEnum('MaterialsRheologyZbar'); Index: /issm/trunk-jpl/test/NightlyRun/test213.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test213.m (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test213.m (revision 16167) @@ -18,5 +18,8 @@ md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy; -md.cluster=generic('name',oshostname(),'np',3); +md.verbose.control=true; +md.verbose.solution=true; + +md.cluster=generic('name',oshostname(),'np',1); md=solve(md,StressbalanceSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test270.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test270.m (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test270.m (revision 16167) @@ -2,8 +2,4 @@ md=setmask(md,'all',''); md=parameterize(md,'../Par/SquareShelf.par'); -md.materials=matdamageice(); -md.materials.rheology_B=paterson(md.initialization.temperature); -md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); -md.damage.isdamage=1; md.damage.D=0.5*ones(md.mesh.numberofvertices,1); md=setflowequation(md,'SSA','all'); Index: /issm/trunk-jpl/test/NightlyRun/test270.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test270.py (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test270.py (revision 16167) @@ -4,5 +4,4 @@ from setmask import * from parameterize import * -from matdamageice import * from paterson import * from setflowequation import * @@ -14,8 +13,5 @@ md=setmask(md,'all','') md=parameterize(md,'../Par/SquareShelf.py') -md.materials=matdamageice() -md.materials.rheology_B=paterson(md.initialization.temperature) -md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) -md.materials.damage_D=0.5*numpy.ones((md.mesh.numberofvertices,1)) +md.damage.D=0.5*numpy.ones((md.mesh.numberofvertices,1)) md=setflowequation(md,'SSA','all') md.cluster=generic('name',oshostname(),'np',3) Index: /issm/trunk-jpl/test/NightlyRun/test272.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test272.m (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test272.m (revision 16167) @@ -2,8 +2,4 @@ md=setmask(md,'all',''); md=parameterize(md,'../Par/SquareShelf.par'); -md.materials=matdamageice(); -md.materials.rheology_B=paterson(md.initialization.temperature); -md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); -md.damage.isdamage=1; md.damage.D=0.5*ones(md.mesh.numberofvertices,1); md=setflowequation(md,'SSA','all'); @@ -12,5 +8,5 @@ %control parameters md.inversion.iscontrol=1; -md.inversion.control_parameters={'MaterialsDamageDbar'}; +md.inversion.control_parameters={'DamageDbar'}; md.inversion.min_parameters=zeros(md.mesh.numberofvertices,1); md.inversion.max_parameters=(1-10^-13)*ones(md.mesh.numberofvertices,1); @@ -28,10 +24,10 @@ %Fields and tolerances to track changes -field_names ={'Gradient','Misfits','MaterialsDamageDbar','Pressure','Vel','Vx','Vy'}; +field_names ={'Gradient','Misfits','DamageDbar','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.StressbalanceSolution.Gradient1),... (md.results.StressbalanceSolution.J),... - (md.results.StressbalanceSolution.MaterialsDamageDbar),... + (md.results.StressbalanceSolution.DamageDbar),... (md.results.StressbalanceSolution.Pressure),... (md.results.StressbalanceSolution.Vel),... Index: /issm/trunk-jpl/test/NightlyRun/test272.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test272.py (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test272.py (revision 16167) @@ -4,5 +4,4 @@ from setmask import * from parameterize import * -from matdamageice import * from paterson import * from setflowequation import * @@ -14,13 +13,10 @@ md=setmask(md,'all','') md=parameterize(md,'../Par/SquareShelf.py') -md.materials=matdamageice() -md.materials.rheology_B=paterson(md.initialization.temperature) -md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) -md.materials.damage_D=0.5*numpy.ones((md.mesh.numberofvertices,1)) +md.damage.D=0.5*numpy.ones((md.mesh.numberofvertices,1)) md=setflowequation(md,'SSA','all') #control parameters md.inversion.iscontrol=1 -md.inversion.control_parameters=['MaterialsDamageDbar'] +md.inversion.control_parameters=['DamageDbar'] md.inversion.min_parameters=10**-13*numpy.ones((md.mesh.numberofvertices,1)) md.inversion.max_parameters=numpy.ones((md.mesh.numberofvertices,1)) @@ -38,10 +34,10 @@ #Fields and tolerances to track changes -field_names =['Gradient','Misfits','MaterialsRheologyDbar','Pressure','Vel','Vx','Vy'] +field_names =['Gradient','Misfits','DamageDbar','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.StressbalanceSolution.Gradient1,\ md.results.StressbalanceSolution.J,\ - md.results.StressbalanceSolution.MaterialsDamageDbar,\ + md.results.StressbalanceSolution.DamageDbar,\ md.results.StressbalanceSolution.Pressure,\ md.results.StressbalanceSolution.Vel,\ Index: /issm/trunk-jpl/test/NightlyRun/test274.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test274.m (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test274.m (revision 16167) @@ -3,8 +3,4 @@ md=setmask(md,'all',''); md=parameterize(md,'../Par/SquareShelf2.par'); -md.materials=matdamageice(); -md.materials.rheology_B=paterson(md.initialization.temperature); -md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); -md.damage.isdamage=1; md.damage.D=0.5*ones(md.mesh.numberofvertices,1); md=setflowequation(md,'SSA','all'); Index: /issm/trunk-jpl/test/NightlyRun/test274.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test274.py (revision 16166) +++ /issm/trunk-jpl/test/NightlyRun/test274.py (revision 16167) @@ -5,5 +5,4 @@ from setmask import * from parameterize import * -from matdamageice import * from paterson import * from setflowequation import * @@ -16,8 +15,5 @@ md=setmask(md,'all','') md=parameterize(md,'../Par/SquareShelf2.py') -md.materials=matdamageice() -md.materials.rheology_B=paterson(md.initialization.temperature) -md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) -md.materials.damage_D=0.5*numpy.ones((md.mesh.numberofvertices,1)) +md.damage.D=0.5*numpy.ones((md.mesh.numberofvertices,1)) md=setflowequation(md,'SSA','all') Index: /issm/trunk-jpl/test/Par/79North.par =================================================================== --- /issm/trunk-jpl/test/Par/79North.par (revision 16166) +++ /issm/trunk-jpl/test/Par/79North.par (revision 16167) @@ -20,4 +20,7 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); %Friction Index: /issm/trunk-jpl/test/Par/79North.py =================================================================== --- /issm/trunk-jpl/test/Par/79North.py (revision 16166) +++ /issm/trunk-jpl/test/Par/79North.py (revision 16167) @@ -33,4 +33,7 @@ md.initialization.temperature=md.initialization.temperature +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Friction md.friction.coefficient=50.*ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/GiaBenchmarksAB.par =================================================================== --- /issm/trunk-jpl/test/Par/GiaBenchmarksAB.par (revision 16166) +++ /issm/trunk-jpl/test/Par/GiaBenchmarksAB.par (revision 16167) @@ -44,4 +44,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Friction md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/GiaBenchmarksCD.par =================================================================== --- /issm/trunk-jpl/test/Par/GiaBenchmarksCD.par (revision 16166) +++ /issm/trunk-jpl/test/Par/GiaBenchmarksCD.par (revision 16167) @@ -43,4 +43,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Friction md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/ISMIPA.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPA.par (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPA.par (revision 16167) @@ -16,4 +16,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' boundary conditions for stressbalance model'); %Create node on boundary first (because we cannot use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPA.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPA.py (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPA.py (revision 16167) @@ -19,4 +19,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " boundary conditions for stressbalance model" #Create node on boundary first (because we cannot use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPB.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPB.par (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPB.par (revision 16167) @@ -16,4 +16,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' boundary conditions for stressbalance model'); %Create node on boundary first (because we cannot use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPB.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPB.py (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPB.py (revision 16167) @@ -19,4 +19,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " boundary conditions for stressbalance model" #Create node on boundary first (because we cannot use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPC.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPC.par (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPC.par (revision 16167) @@ -17,4 +17,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' boundary conditions for stressbalance model:'); %Create node on boundary first (because we can not use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPC.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPC.py (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPC.py (revision 16167) @@ -20,4 +20,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " boundary conditions for stressbalance model:" #Create node on boundary first (because we can not use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPD.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPD.par (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPD.par (revision 16167) @@ -16,4 +16,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' boundary conditions for stressbalance model:'); %Create node on boundary first (because we can not use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPD.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPD.py (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPD.py (revision 16167) @@ -19,4 +19,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " boundary conditions for stressbalance model:" #Create node on boundary first (because we can not use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPE.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPE.par (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPE.par (revision 16167) @@ -26,4 +26,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' boundary conditions for stressbalance model:'); %Create node on boundary first (because we can not use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPE.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPE.py (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPE.py (revision 16167) @@ -31,4 +31,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " boundary conditions for stressbalance model:" #Create node on boundary first (because we can not use mesh) Index: /issm/trunk-jpl/test/Par/ISMIPF.par =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPF.par (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPF.par (revision 16167) @@ -17,4 +17,7 @@ md.materials.rheology_n=1.*ones(md.mesh.numberofelements,1); md.materials.rheology_law='None'; + +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); disp(' boundary conditions for stressbalance model'); Index: /issm/trunk-jpl/test/Par/ISMIPF.py =================================================================== --- /issm/trunk-jpl/test/Par/ISMIPF.py (revision 16166) +++ /issm/trunk-jpl/test/Par/ISMIPF.py (revision 16167) @@ -20,4 +20,7 @@ md.materials.rheology_n=1.*numpy.ones((md.mesh.numberofelements,1)) md.materials.rheology_law='None' + +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) print " boundary conditions for stressbalance model" Index: /issm/trunk-jpl/test/Par/Pig.par =================================================================== --- /issm/trunk-jpl/test/Par/Pig.par (revision 16166) +++ /issm/trunk-jpl/test/Par/Pig.par (revision 16167) @@ -26,4 +26,7 @@ md.initialization.temperature=md.initialization.temperature; +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Friction md.friction.coefficient=50*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/Pig.py =================================================================== --- /issm/trunk-jpl/test/Par/Pig.py (revision 16166) +++ /issm/trunk-jpl/test/Par/Pig.py (revision 16167) @@ -37,4 +37,7 @@ md.initialization.temperature=md.initialization.temperature +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Friction md.friction.coefficient=50.*ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par (revision 16166) +++ /issm/trunk-jpl/test/Par/RoundSheetEISMINT.par (revision 16167) @@ -20,4 +20,7 @@ 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); + +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); disp(' creating surface mass balance'); Index: /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py (revision 16166) +++ /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py (revision 16167) @@ -23,4 +23,7 @@ md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices,1)) #to have the same B as the analytical solution md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) + +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) print " creating surface mass balance" Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 16166) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 16167) @@ -61,4 +61,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Surface mass balance and basal melting md.surfaceforcings.mass_balance=-10.*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.py (revision 16166) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.py (revision 16167) @@ -68,4 +68,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Surface mass balance and basal melting md.surfaceforcings.mass_balance=-10.*numpy.ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par (revision 16166) +++ /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par (revision 16167) @@ -24,4 +24,7 @@ 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); + +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); disp(' creating surface mass balance'); Index: /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py (revision 16166) +++ /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py (revision 16167) @@ -28,4 +28,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " creating surface mass balance" smb_max=0.5 #m/yr Index: /issm/trunk-jpl/test/Par/SquareEISMINT.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareEISMINT.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareEISMINT.par (revision 16167) @@ -26,4 +26,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' creating surface mass balance'); md.surfaceforcings.mass_balance=0.2*ones(md.mesh.numberofvertices,1); %0m/a Index: /issm/trunk-jpl/test/Par/SquareEISMINT.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareEISMINT.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareEISMINT.py (revision 16167) @@ -29,4 +29,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " creating surface mass balance" md.surfaceforcings.mass_balance=0.2*numpy.ones((md.mesh.numberofvertices,1)) #0m/a Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.par (revision 16167) @@ -29,4 +29,7 @@ md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); + +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); %Friction Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 16167) @@ -40,4 +40,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Friction md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetShelf.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareSheetShelf.par (revision 16167) @@ -32,4 +32,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Accumulation and melting md.surfaceforcings.mass_balance=10.*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetShelf.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareSheetShelf.py (revision 16167) @@ -43,4 +43,7 @@ md.materials.rheology_n=3.*ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Accumulation and melting md.surfaceforcings.mass_balance=10.*numpy.ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/SquareShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareShelf.par (revision 16167) @@ -29,4 +29,7 @@ md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Friction md.friction.coefficient=20*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/SquareShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareShelf.py (revision 16167) @@ -59,4 +59,7 @@ md.materials.rheology_n = 3.*ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Friction md.friction.coefficient = 20.*ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/SquareShelf2.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf2.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareShelf2.par (revision 16167) @@ -29,4 +29,7 @@ md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Friction md.friction.coefficient=20*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/SquareShelf2.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf2.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareShelf2.py (revision 16167) @@ -59,4 +59,7 @@ md.materials.rheology_n = 3.*ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Friction md.friction.coefficient = 20.*ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelfConstrained.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.par (revision 16167) @@ -29,4 +29,7 @@ md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + %Surface mass balance and basal melting md.surfaceforcings.mass_balance=10*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py (revision 16167) @@ -40,4 +40,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + #Surface mass balance and basal melting md.surfaceforcings.mass_balance=10.*numpy.ones((md.mesh.numberofvertices,1)) Index: /issm/trunk-jpl/test/Par/SquareThermal.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareThermal.par (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareThermal.par (revision 16167) @@ -27,4 +27,7 @@ md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); +%Damage +md.damage.D=zeros(md.mesh.numberofvertices,1); + disp(' creating surface mass balance'); md.surfaceforcings.mass_balance=ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a Index: /issm/trunk-jpl/test/Par/SquareThermal.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareThermal.py (revision 16166) +++ /issm/trunk-jpl/test/Par/SquareThermal.py (revision 16167) @@ -31,4 +31,7 @@ md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) +#Damage +md.damage.D=zeros((md.mesh.numberofvertices,1)) + print " creating surface mass balance" md.surfaceforcings.mass_balance=numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts #1m/a