Changeset 24921 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

Timestamp:

05/29/20 16:03:44 (5 years ago)

Author:

Eric.Larour

Message:

CHG: 10% increase in sea level solver speed, using a G_rigid kernel pre computed
in the SealevelriseAnalysis section, before G_elastic kernels.

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -5598 +5598 @@
 
         /*elastic green function:*/
+        int         index;
         IssmDouble* G_elastic_precomputed=NULL;
+        IssmDouble* G_rigid_precomputed=NULL;
         int         M;
 
@@ -5649 +5651 @@
 
         /*recover elastic green function:*/
+        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
+
         if(computeelastic){
-                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum));
-                _assert_(parameter);
+                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
                 parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
         }
@@ -5748 +5752 @@
                 for(int i=0;i<gsize;i++){
 
-                        IssmDouble G_rigid=0;  //do not remove =0!
                         IssmDouble G_elastic=0;  //do not remove =0!
 
-                        /*Compute alpha angle between centroid and current vertex : */
+                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
                         lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
-
                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-
-                        //Rigid earth gravitational perturbation:
-                        if(computerigid)G_rigid=1.0/2.0/sin(alpha/2.0);
+                        index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
 
                         //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
-                        if(computeelastic){
-                                int index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
-                                G_elastic += G_elastic_precomputed[index];
-                        }
+                        if(computeelastic) G_elastic += G_elastic_precomputed[index];
 
                         /*Add all components to the Sgi or Sgo solution vectors:*/
-                        Sgi[i]+=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
-
+                        Sgi[i]+=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid_precomputed[index]+G_elastic);
                 }
         }
@@ -5940 +5936 @@
 
         /*precomputed elastic green functions:*/
+        IssmDouble* G_rigid_precomputed = NULL;
         IssmDouble* G_elastic_precomputed = NULL;
         int         M;
-
-        /*computation of Green functions:*/
-        IssmDouble* G_elastic= NULL;
-        IssmDouble* G_rigid= NULL;
+        int         index;
+        IssmDouble alpha;
+        IssmDouble delPhi,delLambda;
 
         /*optimization:*/
@@ -6020 +6016 @@
         /*}}}*/
 
+        /*recover rigid and elastic green functions:*/
+        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
+
         if(computeelastic){
-
                 /*recover elastic green function:*/
-                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
+                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
                 parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
-
-                /*initialize G_elastic:*/
-                G_elastic=xNewZeroInit<IssmDouble>(gsize);
-        }
-        if(computerigid) G_rigid=xNewZeroInit<IssmDouble>(gsize);
-
-        IssmDouble alpha;
-        IssmDouble delPhi,delLambda;
-
+        }
+
+        
         for(int i=0;i<gsize;i++){
 
                 /*Compute alpha angle between centroid and current vertex : */
                 lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
-
                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
                 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
+                index=reCast<int,IssmDouble>(alpha/PI*(M-1));
 
                 /*Rigid earth gravitational perturbation: */
                 if(computerigid){
-                        G_rigid[i]=1.0/2.0/sin(alpha/2.0);
-                        //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
-                        Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
+                        Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid_precomputed[index];
                 }
 
                 /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
                 if(computeelastic){
-                        int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
-                        G_elastic[i] += G_elastic_precomputed[index];
-                        //values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
-                        Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
-                }
-        }
-
-
-        /*Free ressources:*/
-        if(computeelastic) xDelete<IssmDouble>(G_elastic);
-        if(computerigid) xDelete<IssmDouble>(G_rigid);
+                        Sgo[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic_precomputed[index];
+                }
+        }
+
 
         return;