Changeset 26119

Timestamp:

03/19/21 12:37:01 (4 years ago)

Author:

Eric.Larour

Message:

CHG: rearranging sea level optimized core.

Location:

issm/trunk-jpl/src/c/cores

Files:

• cores.h (modified) (1 diff)
• sealevelchange_core.cpp (modified) (10 diffs)

issm/trunk-jpl/src/c/cores/cores.h

@@ -63 +63 @@
 #endif
 void grd_core(FemModel* femmodel);
+void grd_core_optim(FemModel* femmodel);
 void solidearthexternal_core(FemModel* femmodel);
 void dynstr_core(FemModel* femmodel);
-SealevelMasks* sealevel_masks(FemModel* femmodel);
 Vector<IssmDouble>* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea);
 Vector<IssmDouble>* sealevelchange_core_sal(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_barystatic,IssmDouble oceanarea);

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

@@ -17 +17 @@
 void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
+SealevelMasks* sealevel_masks(FemModel* femmodel);
 /*}}}*/
 
@@ -390 +391 @@
 }; /*}}}*/
 
-void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
+void grd_core_optim(FemModel* femmodel) { /*{{{*/
 
         /*variables:{{{*/
         int nel;
         BarystaticContributions* barycontrib=NULL;
+        SealevelMasks* masks=NULL;
         GenericParam<BarystaticContributions*>* barycontribparam=NULL;
-        IssmDouble               barystatic;
         
         Vector<IssmDouble>*    loads=NULL;
@@ -402 +403 @@
         Vector<IssmDouble>*    sealevelloads=NULL;
         Vector<IssmDouble>*    oldsealevelloads=NULL;
-        IssmDouble             sealevelloadsaverage;
         IssmDouble*            allsealevelloads=NULL;
         Vector<IssmDouble>*    oceanareas=NULL;
@@ -413 +413 @@
         int                  step;
         IssmDouble           time; 
-        
-        IssmDouble cumbslc;
-        IssmDouble cumbslcice;
-        IssmDouble cumbslchydro;
+
+        int  modelid,earthid;
+        int  horiz;
+        int  count,frequency,iscoupling;
+        int  grd=0;
         /*}}}*/
 
-        /*retrieve parameters: */
+        /*Verbose: */
+        if(VerboseSolution()) _printf0_("         computing GRD patterns\n");
+
+        /*retrieve parameters:{{{*/
+        femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 
+        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
+        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
+        /*}}}*/
+
+        /*only run if grd was requested, if we are the earth, and we have reached
+         * the necessary number of time steps dictated by :*/
+        if(!grd)            return;
+        if(count!=frequency)return;
+        femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
+        if(iscoupling){
+                femmodel->parameters->FindParam(&modelid,ModelIdEnum);
+                femmodel->parameters->FindParam(&earthid,EarthIdEnum);
+                if(modelid!=earthid)return;
+        }
+        /*retrieve parameters: {{{*/ 
         femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
         barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
@@ -426 +446 @@
         femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
         femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
-        
+        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
+        /*}}}*/
+
         /*initialize matrices and vectors:*/
         femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
@@ -432 +454 @@
         sealevelloads=new Vector<IssmDouble>(nel);
         oceanareas=new Vector<IssmDouble>(nel);
+
+        /*call masks core: */
+        masks=sealevel_masks(femmodel);
         
         /*buildup loads: */
@@ -439 +464 @@
         }
 
-        //Communicate loads from local to global: 
+        //broadcast loads 
         loads->Assemble(); allloads=loads->ToMPISerial();
 
@@ -447 +472 @@
                 element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
         }
-
+        sealevelloads->Assemble();
+        
         //Get ocean area: 
         oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
 
-        //Get sea level loads ocean average: 
-        sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
-        
         //substract ocean average and barystatic contributionfrom sea level loads: 
-        barystatic=barycontrib->Total()/oceanarea/rho_water;
-        sealevelloads->Shift(-sealevelloadsaverage+barystatic);
+        sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
+        
+        //broadcast sea level loads 
         allsealevelloads=sealevelloads->ToMPISerial();
 
@@ -471 +495 @@
                 }
                 sealevelloads->Assemble();
+        
                 
                 //substract ocean average and barystatic contribution
-                sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
-                sealevelloads->Shift(-sealevelloadsaverage+barystatic);
-
-                //broadcast loads 
+                sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
+                
+                //broadcast sea level loads 
                 allsealevelloads=sealevelloads->ToMPISerial();
 
@@ -484 +508 @@
         }
 
+        /*convolve loads and sea level loads to get the deformation:*/
+        for(Object* & object : femmodel->elements->objects){
+                Element* element = xDynamicCast<Element*>(object);
+                //element->SealevelchangeDeformationConvolution(allsealevelloads, allloads,masks);
+        }
+
+        /*Update bedrock motion and geoid:*/
+        /*inputs->AXPY(SealevelEnum,SealevelGRD,1);
+        inputs->AXPY(BedEnum,BedGRD,1);
+        if(horiz){
+                inputs->AXPY(BedrockeastEnum,BedEastGRD,1);
+                inputs->AXPY(BedrocknorthEnum,BedNorthGRD,1);
+        }*/
+
         /*cumulate barystatic contributions and save to results: */
         barycontrib->Cumulate(femmodel->parameters);
         barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
+
 
 }