source: issm/oecreview/Archive/25834-26739/ISSM-26118-26119.diff@ 26740

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

@@ -62 +62 @@
 void sealevelchange_geometry(FemModel* femmodel);
 #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);
 void sealevelchange_core_deformation(Vector<IssmDouble>** pN_radial, Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks);

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

@@ -16 +16 @@
 void TransferSealevel(FemModel* femmodel,int forcingenum);
 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);
 /*}}}*/
 
 /*main cores:*/
@@ -389 +390 @@
         }
 }; /*}}}*/
 
-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;
         IssmDouble*            allloads=NULL; 
         Vector<IssmDouble>*    sealevelloads=NULL;
         Vector<IssmDouble>*    oldsealevelloads=NULL;
-        IssmDouble             sealevelloadsaverage;
         IssmDouble*            allsealevelloads=NULL;
         Vector<IssmDouble>*    oceanareas=NULL;
         IssmDouble             oceanarea;
@@ -412 +412 @@
         IssmDouble           eps_abs;
         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));
         barycontrib=barycontribparam->GetParameterValue();
@@ -425 +445 @@
         femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
         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);
         loads=new Vector<IssmDouble>(nel);
         sealevelloads=new Vector<IssmDouble>(nel);
         oceanareas=new Vector<IssmDouble>(nel);
+
+        /*call masks core: */
+        masks=sealevel_masks(femmodel);
         
         /*buildup loads: */
         for(Object* & object : femmodel->elements->objects){
@@ -438 +463 @@
                 element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); 
         }
 
-        //Communicate loads from local to global: 
+        //broadcast loads 
         loads->Assemble(); allloads=loads->ToMPISerial();
 
         /*convolve loads:*/
@@ -446 +471 @@
                 Element* element = xDynamicCast<Element*>(object);
                 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: 
+        sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
         
-        //substract ocean average and barystatic contributionfrom sea level loads: 
-        barystatic=barycontrib->Total()/oceanarea/rho_water;
-        sealevelloads->Shift(-sealevelloadsaverage+barystatic);
+        //broadcast sea level loads 
         allsealevelloads=sealevelloads->ToMPISerial();
 
         bool converged=false;
@@ -470 +494 @@
                         element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
                 }
                 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();
 
                 //convergence?
@@ -483 +507 @@
                 if (converged)break;
         }
 
+        /*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);
 
+
 }
 /*}}}*/