Changeset 26075

Timestamp:

03/11/21 10:27:18 (4 years ago)

Author:

Eric.Larour

Message:

CHG: removed GIA rates from the test2001 archives.
Transferred new SolidEarthUpdates logic to tws mass transport.

Location:

issm/trunk-jpl

Files:

• src/c/analyses/HydrologyTwsAnalysis.cpp (modified) (3 diffs)
• src/c/cores/hydrology_core.cpp (modified) (4 diffs)
• src/c/cores/masstransport_core.cpp (modified) (3 diffs)
• test/NightlyRun/test2001.m (modified) (1 diff)

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

@@ -38 +38 @@
 
         /*Initialize sea level cumulated sea level loads :*/
-        InputUpdateFromConstantx(inputs,elements,0.,AccumulatedDeltaTwsEnum);
-        InputUpdateFromConstantx(inputs,elements,0.,OldAccumulatedDeltaTwsEnum);
+        iomodel->ConstantToInput(inputs,elements,0.,AccumulatedDeltaTwsEnum,P1Enum);
+        iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaTwsEnum,P1Enum);
         iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum);
 
@@ -79 +79 @@
         if(!element->IsOnBase()) return;
 
-        /*deal with logic of accumulating thickness if we are coupled to a 
-         * sea level core:*/
-        int frequency,count,isgrd;
-        element->FindParam(&isgrd,SolidearthSettingsGRDEnum);
-        if(isgrd){
-                element->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
-                element->FindParam(&count,SealevelchangeRunCountEnum);
-        }       
         /*Fetch dof list and allocate solution vector*/
         int *doflist = NULL;
@@ -106 +98 @@
         xDelete<IssmDouble>(watercolumn);
 
-        /*Get basal element*/
-        int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
-        Element* basalelement=element;
-        if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
-
-        /*Fetch number of nodes and dof for this finite element*/
-        int numvertices = basalelement->GetNumberOfVertices();
-
-        /*Now, we need to do some "processing"*/
-        watercolumn  = xNew<IssmDouble>(numvertices);
-        IssmDouble* oldwatercolumn      = xNew<IssmDouble>(numvertices);
-        IssmDouble* cumdeltawatercolumn = xNew<IssmDouble>(numvertices);
-        IssmDouble* oldcumdeltawatercolumn = xNew<IssmDouble>(numvertices);
-        IssmDouble* deltawatercolumn = xNew<IssmDouble>(numvertices);
-
-        /*Get previous base, tws, surfac and current sealevel and bed:*/
-        basalelement->GetInputListOnVertices(&watercolumn[0],WatercolumnEnum);
-        basalelement->GetInputListOnVertices(&oldwatercolumn[0],WaterColumnOldEnum);
-        basalelement->GetInputListOnVertices(&cumdeltawatercolumn[0],AccumulatedDeltaTwsEnum);
-
-        /*What is the delta tws forcing the sea-level change core: cumulated over time, hence the +=:*/
-        if(isgrd){
-                for(int i=0;i<numvertices;i++){
-                        cumdeltawatercolumn[i] += watercolumn[i]-oldwatercolumn[i];
-                }
-        }
-
-        /*Add input to the element: */
-        if(isgrd){
-                element->AddBasalInput(AccumulatedDeltaTwsEnum,cumdeltawatercolumn,P1Enum);
-                if(count==frequency){
-                        basalelement->GetInputListOnVertices(&oldcumdeltawatercolumn[0],OldAccumulatedDeltaTwsEnum);
-                        element->AddBasalInput(OldAccumulatedDeltaTwsEnum,cumdeltawatercolumn,P1Enum);
-                        for(int i=0;i<numvertices;i++)deltawatercolumn[i]=cumdeltawatercolumn[i]-oldcumdeltawatercolumn[i];
-                        element->AddBasalInput(DeltaTwsEnum,deltawatercolumn,P1Enum);
-                }
-        }
-
-        /*Free ressources:*/
-        xDelete<IssmDouble>(watercolumn);
-        xDelete<IssmDouble>(deltawatercolumn);
-        xDelete<IssmDouble>(cumdeltawatercolumn);
-        xDelete<IssmDouble>(oldcumdeltawatercolumn);
-        xDelete<IssmDouble>(oldwatercolumn);
-        xDelete<int>(doflist);
-        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void           HydrologyTwsAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

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

@@ -9 +9 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
-
-void hydrology_core(FemModel* femmodel){
+void SolidEarthWaterUpdates(FemModel* femmodel);
+
+void hydrology_core(FemModel* femmodel){ /*{{{*/
 
         /*Start profiler*/
@@ -43 +44 @@
                 /*transfer water column thickness to old water column thickness: */
                 InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
+                /*solid earth considerations:*/
+                SolidEarthWaterUpdates(femmodel);
+
 
         }
@@ -59 +63 @@
                 GetVectorFromInputsx(&ug,femmodel,HydrologyTwsSpcEnum,VertexPIdEnum);
                 InputUpdateFromSolutionx(femmodel,ug); 
+
+                /*solid earth considerations:*/
+                SolidEarthWaterUpdates(femmodel);
+
                 delete ug;
         
@@ -219 +227 @@
         /*End profiler*/
         femmodel->profiler->Stop(HYDROLOGYCORE);
-}
+} /*}}}*/
+void SolidEarthWaterUpdates(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, WaterColumnOldEnum,WatercolumnEnum,DeltaTwsEnum);
+        femmodel->inputs->AXPY(+1, DeltaTwsEnum,AccumulatedDeltaTwsEnum,DummyEnum);
+        femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaTwsEnum);
+
+        /*compute total water column change between two sea-level solver time steps, ie. every frequency*dt:*/
+        if(count==frequency){
+                femmodel->inputs->AXPY(-1, OldAccumulatedDeltaTwsEnum,AccumulatedDeltaTwsEnum,DeltaTwsEnum);
+                femmodel->inputs->DuplicateInput(AccumulatedDeltaTwsEnum,OldAccumulatedDeltaTwsEnum);
+        }
+        return;
+}/*}}}*/

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

@@ -10 +10 @@
 #include "../solutionsequences/solutionsequences.h"
 #include "../classes/Inputs/TransientInput.h"
-void SolidEarthUpdates(FemModel* femmodel);
-
+void SolidEarthIceUpdates(FemModel* femmodel);
 void masstransport_core(FemModel* femmodel){ /*{{{*/
 
@@ -77 +76 @@
                         //}
                 }
-                SolidEarthUpdates(femmodel);
+                SolidEarthIceUpdates(femmodel);
                 
                 femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
@@ -100 +99 @@
         femmodel->profiler->Stop(MASSTRANSPORTCORE);
 } /*}}}*/
-void SolidEarthUpdates(FemModel* femmodel){ /*{{{*/
+void SolidEarthIceUpdates(FemModel* femmodel){ /*{{{*/
 
         int isgrd;

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

@@ -50 +50 @@
 
 %Fields and tolerances to track changes
-field_names     ={'UGia','UGiaRate'};
-field_tolerances={1e-13,1e-13};
-field_values={...
-        (md.results.TransientSolution(2).SealevelUGrd),
-        (md.results.TransientSolution(2).SealevelUGrd)
-        };
+field_names     ={'UGrd'};
+field_tolerances={1e-13};
+field_values={md.results.TransientSolution(2).SealevelUGrd};