Changeset 26073

Timestamp:

03/11/21 10:05:44 (4 years ago)

Author:

Eric.Larour

Message:

CHG: fixed test2001 by modifying the mass transport core. Now create inputs
directly in the core, not the MassTransportAnalysis::InputUpdateFromSolution.
Also need a new IoModel::InputToConstant routine.

Location:

issm/trunk-jpl

Files:

• src/c/analyses/MasstransportAnalysis.cpp (modified) (5 diffs)
• src/c/classes/Elements/Element.cpp (modified) (2 diffs)
• src/c/classes/Elements/Element.h (modified) (3 diffs)
• src/c/classes/Elements/Penta.h (modified) (1 diff)
• src/c/classes/Elements/Seg.h (modified) (1 diff)
• src/c/classes/Elements/Tetra.h (modified) (1 diff)
• src/c/classes/Elements/Tria.cpp (modified) (1 diff)
• src/c/classes/Elements/Tria.h (modified) (1 diff)
• src/c/classes/IoModel.cpp (modified) (1 diff)
• src/c/classes/IoModel.h (modified) (1 diff)
• src/c/cores/masstransport_core.cpp (modified) (4 diffs)
• src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp (modified) (2 diffs)
• src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h (modified) (1 diff)
• test/NightlyRun/test2001.m (modified) (3 diffs)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -227 +227 @@
 
         /*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 +860 @@
         newthickness  = xNew<IssmDouble>(numvertices);
         IssmDouble* oldthickness      = xNew<IssmDouble>(numvertices);
-        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
-        IssmDouble* deltathickness = xNew<IssmDouble>(numvertices);
         IssmDouble* newbase           = xNew<IssmDouble>(numvertices);
         IssmDouble* bed               = xNew<IssmDouble>(numvertices);
@@ -877 +875 @@
         basalelement->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
         basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum);
-        if(isgrd)basalelement->GetInputListOnVertices(&cumdeltathickness[0],AccumulatedDeltaIceThicknessEnum);
 
         /*Do we do grounding line migration?*/
@@ -883 +880 @@
         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;i<numvertices;i++){
-                        cumdeltathickness[i] += (newthickness[i]-oldthickness[i]);
-                }
-        }
 
         /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
@@ -924 +914 @@
         element->AddBasalInput(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<int>(numvertices);
-                        vertexsids      = xNew<int>(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<int>(vertexlids);
-                        xDelete<int>(vertexsids);
-                }
-                        
-                /*compute total ice thickness change between two sea-level solver time steps, ie. every frequency*dt:*/
-                if(count==frequency){
-                        IssmDouble* oldcumdeltathickness = xNew<IssmDouble>(numvertices);
-
-                        basalelement->GetInputListOnVertices(&oldcumdeltathickness[0],OldAccumulatedDeltaIceThicknessEnum);
-                        element->AddBasalInput(OldAccumulatedDeltaIceThicknessEnum,cumdeltathickness,P1Enum);
-                        for(int i=0;i<numvertices;i++)deltathickness[i]=cumdeltathickness[i]-oldcumdeltathickness[i];
-                        element->AddBasalInput(DeltaIceThicknessEnum,deltathickness,P1Enum);
-
-                        xDelete<IssmDouble>(oldcumdeltathickness);
-                }
-
-        }
+        
 
         /*Free ressources:*/
         xDelete<IssmDouble>(newthickness);
-        xDelete<IssmDouble>(cumdeltathickness);
-        xDelete<IssmDouble>(deltathickness);
         xDelete<IssmDouble>(newbase);
         xDelete<IssmDouble>(newsurface);

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

@@ -1709 +1709 @@
 }
 /*}}}*/
+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<NUM_VERTICES;i++){
+                vertexids[i] =reCast<int>(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;i<NUM_VERTICES;i++) values[i]=value_in;
+        this->SetElementInput(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 +1966 @@
         /*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);
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -130 +130 @@
                 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 +136 @@
                 void               InputUpdateFromConstant(int constant, int name);
                 void               InputUpdateFromConstant(bool constant, int name);
+                void               InputUpdateFromConstant(IssmDouble constant, int name, int type);
 
                 bool               IsFloating();
@@ -331 +333 @@
                 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");};

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -165 +165 @@
                 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);

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -139 +139 @@
            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");};

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -143 +143 @@
                 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);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -4177 +4177 @@
         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");
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -126 +126 @@
                 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);

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -447 +447 @@
 
         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<Element*>(object);
+                        element->InputCreateP1FromConstant(inputs,this,value,vector_enum);
+                }
+        }
+        else _error_("not supported yet!");
+        return;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/IoModel.h

@@ -123 +123 @@
                 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);

issm/trunk-jpl/src/c/cores/masstransport_core.cpp

@@ -9 +9 @@
 #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 +77 @@
                         //}
                 }
+                SolidEarthUpdates(femmodel);
+                
                 femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
                 extrudefrombase_core(femmodel);
@@ -81 +85 @@
                 femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
                 extrudefrombase_core(femmodel);
+                
         }
 
@@ -94 +99 @@
         /*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<Element*>(object);
+                        const int numvertices = element->GetNumberOfVertices();
+                        IssmDouble* cumdeltathickness=NULL;
+
+                        /*Get values and lid list and recover vertices ids needed to initialize inputs*/
+                        vertexlids      = xNew<int>(numvertices);
+                        vertexsids      = xNew<int>(numvertices);
+                        cumdeltathickness=xNew<IssmDouble>(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<int>(vertexlids);
+                        xDelete<int>(vertexsids);
+                        xDelete<IssmDouble>(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;
+}/*}}}*/

issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp

@@ -27 +27 @@
         }
 }
+
+
+
 void InputUpdateFromConstantx(FemModel* femmodel,IssmDouble constant, int name){
 
@@ -60 +63 @@
         }
 }
+
+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<Element*>(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){
 

issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h

@@ -17 +17 @@
 #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);
 

issm/trunk-jpl/test/NightlyRun/test2001.m

@@ -7 +7 @@
 %GIA Ivins, 2 layer model.
 md.solidearth.settings.grdmodel=2;
+md.solidearth.settings.isgrd=1;
 
 md.materials=materials('litho','ice');
@@ -22 +23 @@
 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 +50 @@
 
 %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)
         };