Changeset 24916

Timestamp:

05/28/20 23:43:54 (5 years ago)

Author:

Eric.Larour

Message:

CHG: speeding up by not recomputing areas all the time.

Location:

issm/trunk-jpl/src/c

Files:

• classes/FemModel.cpp (modified) (10 diffs)
• classes/FemModel.h (modified) (1 diff)
• cores/cores.h (modified) (1 diff)
• cores/sealevelrise_core.cpp (modified) (11 diffs)

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

@@ -4684 +4684 @@
 #endif
 #ifdef _HAVE_SEALEVELRISE_
-void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
+void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
 
         /*serialized vectors:*/
@@ -4706 +4706 @@
         /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
         for(int i=0;i<elements->Size();i++){
+                IssmDouble area;
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                oceanarea_cpu += element->OceanArea();
-                eartharea_cpu += element->GetAreaSpherical();
+                area=element->GetAreaSpherical();
+                eartharea_cpu += area;
+                if (element->IsOceanInElement()) oceanarea_cpu += area;
         }
         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
@@ -4737 +4739 @@
 
         /*Assign output pointers:*/
+        *peartharea=eartharea;
+        *poceanarea=oceanarea;
         *peustatic=eustatic;
 
 }
 /*}}}*/
-void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
+void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
 
         /*serialized vectors:*/
         IssmDouble* RSLg_old=NULL;
-
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
 
         IssmDouble* RSLgo  = NULL;
@@ -4762 +4763 @@
         RSLg_old=pRSLg_old->ToMPISerial();
 
-        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
-        for(int i=0;i<elements->Size();i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
         /*Call the sea level rise core: */
         for(int i=0;i<elements->Size();i++){
@@ -4786 +4778 @@
 }
 /*}}}*/
-void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
 
         /*serialized vectors:*/
         IssmDouble* RSLg_old=NULL;
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
         IssmDouble      tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g;
         IssmDouble      load_love_k2 = -0.30922675; //degree 2 load Love number
@@ -4799 +4789 @@
         /*Serialize vectors from previous iteration:*/
         RSLg_old=pRSLg_old->ToMPISerial();
-
-        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
-        for(int i=0;i<elements->Size();i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
         IssmDouble moi_list[3]={0,0,0};
@@ -4872 +4854 @@
 }
 /*}}}*/
-void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
+void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
 
         /*serialized vectors:*/
         IssmDouble* RSLg=NULL;
-
-        IssmDouble  eartharea=0;
-        IssmDouble  eartharea_cpu=0;
 
         IssmDouble* Up  = NULL;
@@ -4898 +4877 @@
         indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
 
-        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
-        for(int i=0;i<elements->Size();i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                eartharea_cpu += element->GetAreaSpherical();
-        }
-        ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
         /*Call the sea level rise core: */
         for(int i=0;i<elements->Size();i++){
@@ -4929 +4900 @@
 }
 /*}}}*/
-IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg) { /*{{{*/
+IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg,IssmDouble oceanarea) { /*{{{*/
 
         IssmDouble* RSLg_serial=NULL;
         IssmDouble  oceanvalue,oceanvalue_cpu;
-        IssmDouble  oceanarea,oceanarea_cpu;
 
         /*Serialize vectors from previous iteration:*/
@@ -4940 +4910 @@
         /*Initialize:*/
         oceanvalue_cpu=0;
-        oceanarea_cpu=0;
 
         /*Go through elements, and add contribution from each element and divide by overall ocean area:*/
         for(int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                oceanarea_cpu += element->OceanArea();
                 oceanvalue_cpu += element->OceanAverage(RSLg_serial);
         }
-        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-        //_printf0_("Ocean area: " << oceanarea << "\n");
-
+        
         ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

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

@@ -162 +162 @@
                 #endif
                 #ifdef _HAVE_SEALEVELRISE_
-                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
-                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
-                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
-                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
-                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
+                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
+                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
+                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
+                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
+                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg,IssmDouble oceanarea);
                 #endif
                 void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);

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

@@ -55 +55 @@
 void geodetic_core(FemModel* femmodel);
 void steric_core(FemModel* femmodel);
-Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel);
-Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic);
-void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg);
+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea);
+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea);
+void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea);
 void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia,Vector<IssmDouble>** pN_gia,FemModel*  femmodel,Vector<IssmDouble>* RSLg);
 void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg);

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

@@ -92 +92 @@
         int  geodetic=0;
         IssmDouble          dt;
+        IssmDouble          eartharea,oceanarea;
 
         /*Should we even be here?:*/
@@ -145 +146 @@
 
                 /*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
-                RSLg_eustatic=sealevelrise_core_eustatic(femmodel); 
+                RSLg_eustatic=sealevelrise_core_eustatic(femmodel,&eartharea,&oceanarea); 
 
                 /*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
-                RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic); 
+                RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic,eartharea,oceanarea); 
 
                 /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
-                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg);
+                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea);
 
                 /*compute viscosus (GIA) geodetic signatures:*/
@@ -303 +304 @@
 }
 /*}}}*/
-Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/
+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/
 
         /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
@@ -317 +318 @@
         IssmDouble *radius    = NULL;
         int         loop;
+        IssmDouble eartharea,oceanarea;
 
         /*outputs:*/
@@ -334 +336 @@
 
         /*call the eustatic main module: */
-        femmodel->SealevelriseEustatic(RSLgi,&eustatic, latitude, longitude, radius,loop); //this computes 
+        femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, latitude, longitude, radius,loop); //this computes 
 
         /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
-        RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi);
+        RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,oceanarea);
 
         /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the 
@@ -350 +352 @@
         xDelete<IssmDouble>(longitude);
         xDelete<IssmDouble>(radius);
+
+        /*Assign output pointers and return: */
+        *peartharea=eartharea;
+        *poceanarea=oceanarea;
         return RSLgi;
 }/*}}}*/
-Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic){ /*{{{*/
+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/
 
         /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
@@ -419 +425 @@
 
                 /*call the non eustatic module: */
-                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop);
+                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  eartharea,latitude,longitude, radius,verboseconvolution,loop);
 
                 /*assemble solution vector: */
@@ -428 +434 @@
                         /*call rotational feedback  module: */
                         RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
-                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius); 
+                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,eartharea,latitude,longitude,radius); 
                         RSLgo_rot->Assemble(); 
 
@@ -440 +446 @@
 
                 /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
-                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo);
+                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,oceanarea);
 
                 /*RSLg is the sum of the eustatic term, and the ocean terms: */
@@ -477 +483 @@
         return RSLg;
 } /*}}}*/
-void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
+void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea){ /*{{{*/
 
         Vector<IssmDouble> *U_esa  = NULL; 
@@ -515 +521 @@
 
         /*call the elastic main modlule:*/ 
-        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
+        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea,latitude,longitude,radius,xx,yy,zz,loop,horiz);
 
         /*Assign output pointers:*/