Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26072) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26073) @@ -227,6 +227,6 @@ /*Initialize sea level cumulated sea level loads :*/ - InputUpdateFromConstantx(inputs,elements,0.,AccumulatedDeltaIceThicknessEnum); - InputUpdateFromConstantx(inputs,elements,0.,OldAccumulatedDeltaIceThicknessEnum); + iomodel->ConstantToInput(inputs,elements,0,AccumulatedDeltaIceThicknessEnum,P1Enum); + iomodel->ConstantToInput(inputs,elements,0,OldAccumulatedDeltaIceThicknessEnum,P1Enum); /*for Ivins deformation model, initialize history of ice thickness changes:*/ @@ -860,6 +860,4 @@ newthickness = xNew(numvertices); IssmDouble* oldthickness = xNew(numvertices); - IssmDouble* cumdeltathickness = xNew(numvertices); - IssmDouble* deltathickness = xNew(numvertices); IssmDouble* newbase = xNew(numvertices); IssmDouble* bed = xNew(numvertices); @@ -877,5 +875,4 @@ basalelement->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum); basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum); - if(isgrd)basalelement->GetInputListOnVertices(&cumdeltathickness[0],AccumulatedDeltaIceThicknessEnum); /*Do we do grounding line migration?*/ @@ -883,11 +880,4 @@ element->FindParam(&isgroundingline,TransientIsgroundinglineEnum); if(isgroundingline) basalelement->GetInputListOnVertices(&bed[0],BedEnum); - - /*What is the delta thickness forcing the sea-level change core: cumulated over time, hence the +=:*/ - if(isgrd){ - for(int i=0;iAddBasalInput(SurfaceEnum,newsurface,P1Enum); element->AddBasalInput(BaseEnum,newbase,P1Enum); - if(isgrd){ - int grdmodel; - - /*accumulate ice thickness changes: */ - element->AddBasalInput(AccumulatedDeltaIceThicknessEnum,cumdeltathickness,P1Enum); - - /*for Ivins deformation model, keep history of ice thickness changes:*/ - element->FindParam(&grdmodel,GrdModelEnum); - if(grdmodel==IvinsEnum){ - int *vertexlids = NULL; - int *vertexsids= NULL; - IssmDouble time; - - TransientInput* transientinput = basalelement->inputs->GetTransientInput(TransientAccumulatedDeltaIceThicknessEnum); - - /*Get values and lid list and recover vertices ids needed to initialize inputs*/ - vertexlids = xNew(numvertices); - vertexsids = xNew(numvertices); - element->FindParam(&time,TimeEnum); - element->GetVerticesLidList(&vertexlids[0]); - element->GetVerticesSidList(&vertexsids[0]); - - /*Add the current time cumdeltathickness to the existing time series: */ - switch(element->ObjectEnum()){ - case TriaEnum: transientinput->AddTriaTimeInput( time,numvertices,vertexlids,cumdeltathickness,P1Enum); break; - default: _error_("Not implemented yet"); - } - xDelete(vertexlids); - xDelete(vertexsids); - } - - /*compute total ice thickness change between two sea-level solver time steps, ie. every frequency*dt:*/ - if(count==frequency){ - IssmDouble* oldcumdeltathickness = xNew(numvertices); - - basalelement->GetInputListOnVertices(&oldcumdeltathickness[0],OldAccumulatedDeltaIceThicknessEnum); - element->AddBasalInput(OldAccumulatedDeltaIceThicknessEnum,cumdeltathickness,P1Enum); - for(int i=0;iAddBasalInput(DeltaIceThicknessEnum,deltathickness,P1Enum); - - xDelete(oldcumdeltathickness); - } - - } + /*Free ressources:*/ xDelete(newthickness); - xDelete(cumdeltathickness); - xDelete(deltathickness); xDelete(newbase); xDelete(newsurface); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26073) @@ -1709,4 +1709,24 @@ } /*}}}*/ +void Element::InputCreateP1FromConstant(Inputs* inputs,IoModel* iomodel,IssmDouble value_in,int vector_enum){/*{{{*/ + + const int NUM_VERTICES = this->GetNumberOfVertices(); + + int vertexids[MAXVERTICES]; + int vertexlids[MAXVERTICES]; + IssmDouble values[MAXVERTICES]; + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(int i=0;i(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab + vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1]; + } + + for(int i=0;iSetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum); + +} +/*}}}*/ void Element::ControlInputCreate(IssmDouble* vector,IssmDouble* min_vector,IssmDouble* max_vector,Inputs* inputs,IoModel* iomodel,int M,int N,IssmDouble scale,int input_enum,int id){/*{{{*/ @@ -1946,4 +1966,13 @@ /*update input*/ this->SetElementInput(name,constant); +} +/*}}}*/ +void Element::InputUpdateFromConstant(IssmDouble constant, int name, int type){/*{{{*/ + + /*Check that name is an element input*/ + if(!IsInputEnum(name)) return; + + /*update input*/ + this->SetElementInput(name,constant,type); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26073) @@ -130,4 +130,5 @@ int Id(); void InputCreate(IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); + void InputCreateP1FromConstant(Inputs* inputs,IoModel* iomodel,IssmDouble value,int vector_enum); void ControlInputCreate(IssmDouble* doublearray,IssmDouble* independents_min,IssmDouble* independents_max,Inputs*inputs,IoModel* iomodel,int M,int N,IssmDouble scale,int input_enum,int id); void DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum); @@ -135,4 +136,5 @@ void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); + void InputUpdateFromConstant(IssmDouble constant, int name, int type); bool IsFloating(); @@ -331,5 +333,6 @@ virtual void ResetHooks()=0; virtual void RignotMeltParameterization(void){_error_("not implemented yet");}; - virtual void SetElementInput(int enum_in,IssmDouble values){_error_("not implemented yet");}; + virtual void SetElementInput(int enum_in,IssmDouble value){_error_("not implemented yet");}; + virtual void SetElementInput(int enum_in,IssmDouble value, int type)=0; virtual void SetElementInput(Inputs* inputs,int enum_in,IssmDouble values){_error_("not implemented yet");}; virtual void SetElementInput(Inputs* inputs,int numindices,int* indices,IssmDouble* values,int enum_in){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26073) @@ -165,6 +165,8 @@ void ResetHooks(); void RignotMeltParameterization(); - void SetElementInput(int enum_in,IssmDouble values); - void SetElementInput(Inputs* inputs,int enum_in,IssmDouble values); + void SetElementInput(int enum_in,IssmDouble value); + void SetElementInput(int enum_in,IssmDouble value,int type){_error_("not implemented yet");}; + void SetElementInput(Inputs* inputs,int enum_in,IssmDouble value); + void SetElementInput(Inputs* inputs,int enum_in,IssmDouble value,int type){_error_("not implemented yet");}; void SetElementInput(Inputs* inputs,int numindices,int* indices,IssmDouble* values,int enum_in); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int M,int N); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26073) @@ -139,4 +139,5 @@ Element* SpawnBasalElement(bool depthaverage_materials){_error_("not implemented yet");}; Element* SpawnTopElement(void){_error_("not implemented yet");}; + void SetElementInput(int enum_in,IssmDouble value, int type){_error_("not implemented yet");}; IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; void StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26073) @@ -143,4 +143,6 @@ void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");}; + void SetElementInput(Inputs* inputs,int enum_in,IssmDouble value,int type){_error_("not implemented yet");}; + void SetElementInput(int enum_in, IssmDouble value, int type){_error_("not implemented yet");}; Element* SpawnBasalElement(bool depthaverage_materials); Element* SpawnTopElement(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26073) @@ -4177,4 +4177,19 @@ inputs->SetTriaInput(enum_in,P0Enum,this->lid,value); +} +/*}}}*/ +void Tria::SetElementInput(int enum_in,IssmDouble value,int type){/*{{{*/ + + if(type==P0Enum){ + this->inputs->SetTriaInput(enum_in,P0Enum,this->lid,value); + } + else if(type==P1Enum){ + IssmDouble values[3]; + for(int i=0;i<3;i++)values[i]=value; + int lidlist[3]; + this->GetVerticesLidList(&lidlist[0]); + this->inputs->SetTriaInput(enum_in,P1Enum,3,&lidlist[0],&values[0]); + } + else _error_("interpolation type not supported yet"); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26072) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26073) @@ -126,6 +126,7 @@ void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); void RignotMeltParameterization(); - void SetElementInput(int enum_in,IssmDouble values); - void SetElementInput(Inputs* inputs,int enum_in,IssmDouble values); + void SetElementInput(int enum_in,IssmDouble value); + void SetElementInput(int enum_in,IssmDouble value,int type); + void SetElementInput(Inputs* inputs,int enum_in,IssmDouble value); void SetElementInput(Inputs* inputs,int numindices,int* indices,IssmDouble* values,int enum_in); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int M,int N); Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26072) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26073) @@ -447,4 +447,16 @@ return NULL; +} +/*}}}*/ +void IoModel::ConstantToInput(Inputs* inputs,Elements* elements,IssmDouble value, int vector_enum,int type){/*{{{*/ + + if (type==P1Enum){ + for(Object* & object : elements->objects){ + Element* element=xDynamicCast(object); + element->InputCreateP1FromConstant(inputs,this,value,vector_enum); + } + } + else _error_("not supported yet!"); + return; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 26072) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 26073) @@ -123,4 +123,5 @@ void CheckFile(void); Param *CopyConstantObject(const char* constant_name,int param_enum); + void ConstantToInput(Inputs* inputs,Elements* elements,IssmDouble value, int vector_enum,int type); IssmDouble *Data(const char* data_name); void DeclareIndependents(bool trace,IssmPDouble* X); Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26072) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26073) @@ -9,6 +9,8 @@ #include "../modules/modules.h" #include "../solutionsequences/solutionsequences.h" +#include "../classes/Inputs/TransientInput.h" +void SolidEarthUpdates(FemModel* femmodel); -void masstransport_core(FemModel* femmodel){ +void masstransport_core(FemModel* femmodel){ /*{{{*/ /*Start profiler*/ @@ -75,4 +77,6 @@ //} } + SolidEarthUpdates(femmodel); + femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); @@ -81,4 +85,5 @@ femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); + } @@ -94,3 +99,63 @@ /*profiler*/ femmodel->profiler->Stop(MASSTRANSPORTCORE); -} +} /*}}}*/ +void SolidEarthUpdates(FemModel* femmodel){ /*{{{*/ + + int isgrd; + int grdmodel; + IssmDouble time; + int frequency,count; + + /*retrieve parameters:*/ + femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum); + femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); + femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); + + /*early return?:*/ + if(!isgrd)return; + + /*From old and new thickness, create delta thickness and accumulate:*/ + femmodel->inputs->AXPY(-1, ThicknessOldEnum,ThicknessEnum,DeltaIceThicknessEnum); + femmodel->inputs->AXPY(+1, DeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DummyEnum); + femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaIceThicknessEnum); + + /*for Ivins deformation model, keep history of ice thickness changes inside TransientAccumulatedDeltaIceThicknessEnum:*/ + if(grdmodel==IvinsEnum){ + + TransientInput* transientinput = femmodel->inputs->GetTransientInput(TransientAccumulatedDeltaIceThicknessEnum); + + for(Object* & object : femmodel->elements->objects){ + int *vertexlids = NULL; + int *vertexsids= NULL; + Element* element=xDynamicCast(object); + const int numvertices = element->GetNumberOfVertices(); + IssmDouble* cumdeltathickness=NULL; + + /*Get values and lid list and recover vertices ids needed to initialize inputs*/ + vertexlids = xNew(numvertices); + vertexsids = xNew(numvertices); + cumdeltathickness=xNew(numvertices); + element->GetVerticesLidList(&vertexlids[0]); + element->GetVerticesSidList(&vertexsids[0]); + element->GetInputListOnVertices(&cumdeltathickness[0],AccumulatedDeltaIceThicknessEnum); + + /*Add the current time cumdeltathickness to the existing time series: */ + switch(element->ObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( time,numvertices,vertexlids,cumdeltathickness,P1Enum); break; + default: _error_("Not implemented yet"); + } + xDelete(vertexlids); + xDelete(vertexsids); + xDelete(cumdeltathickness); + } + } + + /*compute total ice thickness change between two sea-level solver time steps, ie. every frequency*dt:*/ + if(count==frequency){ + femmodel->inputs->AXPY(-1, OldAccumulatedDeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DeltaIceThicknessEnum); + femmodel->inputs->DuplicateInput(AccumulatedDeltaIceThicknessEnum,OldAccumulatedDeltaIceThicknessEnum); + } + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp (revision 26072) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp (revision 26073) @@ -27,4 +27,7 @@ } } + + + void InputUpdateFromConstantx(FemModel* femmodel,IssmDouble constant, int name){ @@ -60,4 +63,27 @@ } } + +void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,IssmDouble constant, int name,int type){ + + if(type==P0Enum) InputUpdateFromConstantx(inputs, elements, constant,name); + else if(type==P1Enum){ + + if(VerboseModule()) _printf0_(" Input updates from constant (P1 version)\n"); + + /*Elements and loads drive the update: */ + if(IsInputEnum(name)){ + for(Object* & object : elements->objects){ + Element* element = xDynamicCast(object); + element->InputUpdateFromConstant(constant,name,P1Enum); + } + } + else{ + _error_("not supported yet"); + } + } + else _error_("InputUpdateFromConstantx error message: type not supported yet!"); + + +} void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,bool constant, int name){ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h (revision 26072) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h (revision 26073) @@ -17,4 +17,5 @@ #endif void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,IssmDouble constant,int name); +void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,IssmDouble constant,int name, int type); void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,bool constant,int name); Index: /issm/trunk-jpl/test/NightlyRun/test2001.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26072) +++ /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26073) @@ -7,4 +7,5 @@ %GIA Ivins, 2 layer model. md.solidearth.settings.grdmodel=2; +md.solidearth.settings.isgrd=1; md.materials=materials('litho','ice'); @@ -22,6 +23,6 @@ md.timestepping.start_time=-2400000; %4,800 kyr :: EVALUATION TIME md.timestepping.time_step= 2400000; %2,400 kyr :: EVALUATION TIME -% to get rid of default final_time, make sure final_time > start_time -md.timestepping.final_time=2400000; %2,400 kyr +% to get rid of default final_time: make sure final_time>start_time +md.timestepping.final_time=2400000; %2,500 kyr md.masstransport.spcthickness=[... [md.geometry.thickness; 0],... @@ -49,7 +50,8 @@ %Fields and tolerances to track changes -field_names ={'UGrd'}; -field_tolerances={1e-13}; +field_names ={'UGia','UGiaRate'}; +field_tolerances={1e-13,1e-13}; field_values={... - (md.results.TransientSolution(1).SealevelUGrd) + (md.results.TransientSolution(2).SealevelUGrd), + (md.results.TransientSolution(2).SealevelUGrd) };