Changeset 22143

Timestamp:

10/05/17 16:48:16 (8 years ago)

Author:

Eric.Larour

Message:

CHG: better marching through time for slr solution

Location:

issm/branches/trunk-larour-NatGeoScience2016/src

Files:

• c/analyses/MasstransportAnalysis.cpp (modified) (4 diffs)
• c/analyses/SealevelriseAnalysis.cpp (modified) (3 diffs)
• c/cores/cores.h (modified) (1 diff)
• c/cores/sealevelrise_core.cpp (modified) (10 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• m/classes/sealevelmodel.m (modified) (1 diff)

issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/MasstransportAnalysis.cpp

@@ -157 +157 @@
         iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
+        
+        /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on 
+         *elements, otherwise we would have used InputUpdateFromConstantx: */
+        counter=0;
+        for(int i=0;i<iomodel->numberofelements;i++){
+                if(iomodel->my_elements[i]){
+                        Element* element=(Element*)elements->GetObjectByOffset(counter);
+                        element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum);
+                        counter++;
+                }
+        }
 
         if(!issmb){
@@ -650 +661 @@
         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
         IssmDouble* newthickness   = xNew<IssmDouble>(numnodes);
-        IssmDouble* deltathickness = xNew<IssmDouble>(numnodes);
+        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numnodes);
+        IssmDouble* deltathickness    = xNew<IssmDouble>(numnodes);
         IssmDouble* newbase        = xNew<IssmDouble>(numnodes);
         IssmDouble* newsurface     = xNew<IssmDouble>(numnodes);
@@ -675 +687 @@
         basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
         basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum);
-
-        /*What is the delta thickness forcing the sea-level rise core: */
-        for(i=0;i<numnodes;i++) deltathickness[i]=newthickness[i]-oldthickness[i];
+        basalelement->GetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
+
+        /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/
+        for(i=0;i<numnodes;i++){
+                cumdeltathickness[i]+=newthickness[i]-oldthickness[i];
+                deltathickness[i]=newthickness[i]-oldthickness[i];
+        }
 
         /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
@@ -702 +718 @@
         }
 
+
         /*Add input to the element: */
         element->AddBasalInput(ThicknessEnum,newthickness,P1Enum);
-        element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
         element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
         element->AddBasalInput(BaseEnum,newbase,P1Enum);
+        element->AddBasalInput(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum);
+        element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
 
         /*Free ressources:*/

issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp

@@ -39 +39 @@
         iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum);
         iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
-        
+
         /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on 
          *elements, otherwise we would have used InputUpdateFromConstantx: */
@@ -88 +88 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
-        parameters->AddObject(iomodel->CopyConstantObject("md.slr.run_frequency",SealevelriseRunFrequencyEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
@@ -97 +96 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
-        parameters->SetParam(1,SealevelriseRunCountEnum); 
 
         iomodel->FetchData(&elastic,"md.slr.elastic");

issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h

@@ -66 +66 @@
 void TransferSealevel(FemModel* femmodel,int forcingenum);
 void EarthMassTransport(FemModel* femmodel);
-bool SlrAccumulateDeltaThickness(FemModel* femmodel);
 
 //solution configuration

issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp

@@ -19 +19 @@
         Vector<IssmDouble> *U_north   = NULL; 
         Vector<IssmDouble> *U_east    = NULL; 
+        Vector<IssmDouble>* cumdeltathickness=NULL; 
         bool save_results,isslr,iscoupler;
         int configuration_type;
@@ -35 +36 @@
         IssmDouble          *zz     = NULL;
         IssmDouble          dt,steric_rate;
-        int                 frequency;
+        int                 frequency,count;
         int  loop;
         int  horiz;
@@ -45 +46 @@
         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
         femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
+        femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum);
+        femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
+
         
         /*first, recover lat,long and radius vectors from vertices: */
@@ -77 +81 @@
         if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
 
-        /*figure out whether deltathickness was computed on the earth (mass trasnport module should 
+        /*figure out whether deltathickness was computed on the earth (mass transport module should 
          * have taken care of it for ice caps:  */
         if(iscoupler){
@@ -88 +92 @@
         }
 
-        /*transfer deltathickness forcing from ice caps to earth model: */
-        if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum);
+        /*transfer cum deltathickness forcing from ice caps to earth model: */
+        if(iscoupler & (count==frequency)) TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
 
         /*call sea-level rise sub cores:*/
-        if(isslr && !SlrAccumulateDeltaThickness(femmodel)){
+        if(isslr & (count==frequency)){
+                        
+                if(VerboseSolution()) _printf0_("   computing SLR core\n");
+
+                if(iscoupler){
+                        /*we have accumulated thicknesses, dump them in deltathcikness: */
+                        GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+                        InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
+                        delete cumdeltathickness;
+                }
 
                 femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
@@ -154 +167 @@
 
         }
-
+                
+        /*transfer sea-level back to ice caps, reset cum thicknesses to 0 and reset counter: */
+        if(iscoupler){
+                if (count==frequency) {
+                        //transfer sea level back to ice caps:
+                        TransferSealevel(femmodel,SealevelEnum);
+
+                        //reset cumdeltathickness  to 0: 
+                        GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+                        cumdeltathickness->Set(0);
+                        InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+                        delete cumdeltathickness;
+
+                        //reset counter to 1:
+                        femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
+                }
+                else{
+                        //increment counter: 
+                        count++;
+                        femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); 
+                }
+        }
+        
         /*Free ressources: */
         delete longitude; 
@@ -163 +198 @@
         delete zz; 
 
-        /*transfer sea-level back to ice caps: */
-        if(iscoupler)TransferSealevel(femmodel,SealevelEnum);
 } 
 /*}}}*/
@@ -260 +293 @@
                                 /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular 
                                  * ice cap: */
-                                forcingglobal->SetValues(M,index,forcingfromcap,INS_VAL);
+                                forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
                                 xDelete<int>(index);
                         }
@@ -406 +439 @@
         Vector<IssmDouble> *newthickness    = NULL;
         Vector<IssmDouble> *deltathickness    = NULL;
+        Vector<IssmDouble> *cumdeltathickness    = NULL;
         int nv;
         
@@ -430 +464 @@
         InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
 
+        /*add to cumulated delta thickness: */
+        GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+        cumdeltathickness->AXPY(deltathickness,1); 
+        InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+
         /*free ressources:*/
         delete oldthickness;
         delete newthickness;
         delete deltathickness;
+        delete cumdeltathickness;
 
 } /*}}}*/
-bool SlrAccumulateDeltaThickness(FemModel* femmodel){ /*{{{*/
-
-        Vector<IssmDouble>* deltathickness=NULL; 
-        Vector<IssmDouble>* cumdeltathickness=NULL; 
-        int frequency=1;
-        int count;
-        int solution_type;
-        
-        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-        if(solution_type==SealevelriseSolutionEnum)return false ;  //we don't accumulate when we run this solution type!
-
-        femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum);
-        femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
-
-        if (count==frequency){
-                femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
-                
-                /*increment delta thickness load: */
-                GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
-                GetVectorFromInputsx(&deltathickness,femmodel,SealevelriseDeltathicknessEnum,VertexSIdEnum);
-
-                cumdeltathickness->AXPY(deltathickness,1);
-
-                /*update deltathickness with the final cumulated thickness: */
-                InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
-                
-                /*reset cumulative thickness to 0: */
-                cumdeltathickness->Set(0);
-                InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
-                
-                delete deltathickness;
-                delete cumdeltathickness;
-
-                return false;
-        }
-        else{
-                count++;
-                femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); //increment counter
-                
-                /*increment cumulative delta thickness load: */
-                GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
-                GetVectorFromInputsx(&deltathickness,femmodel,SealevelriseDeltathicknessEnum,VertexSIdEnum);
-
-                cumdeltathickness->AXPY(deltathickness,1);
-                InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
-                
-                delete deltathickness;
-                delete cumdeltathickness;
-
-                return true;
-        }
-
-} /*}}}*/

issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -66 +66 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.run_frequency",SealevelriseRunFrequencyEnum));
+        parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1));  
         {/*This is specific to ice...*/
                 parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum));

issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m

@@ -58 +58 @@
                                 end
                         end
+                        
+                        %check that run_frequency is the same everywhere: 
+                        for i=1:length(slm.icecaps),
+                                if slm.icecaps{i}.slr.run_frequency~=slm.earth.slr.run_frequency,
+                                        error(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));
+                                end
+                        end
+
+