Changeset 21261 for issm/trunk-jpl/src/c/classes/FemModel.cpp

Timestamp:

10/12/16 12:16:21 (8 years ago)

Author:

adhikari

Message:

NEW: added new ESA capability to compute 3D elastostatic crustal deformation

File:

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

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

@@ -625 +625 @@
                         break;
                 
+                case EsaSolutionEnum:
+                        analyses_temp[numanalyses++]=EsaAnalysisEnum;
+                        break;
+                
                 case SealevelriseSolutionEnum:
                         analyses_temp[numanalyses++]=SealevelriseAnalysisEnum;
@@ -638 +642 @@
 
                 case TransientSolutionEnum:{
-                        bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslr,isgia;
+                        bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslr,isesa,isgia;
                         iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
                         iomodel->FindConstant(&isFS,"md.flowequation.isFS");
@@ -651 +655 @@
                         iomodel->FindConstant(&issmb,"md.transient.issmb");
                         iomodel->FindConstant(&isslr,"md.transient.isslr");
+                        iomodel->FindConstant(&isesa,"md.transient.isesa");
                         iomodel->FindConstant(&isgia,"md.transient.isgia");
                         if(isstressbalance){
@@ -694 +699 @@
                         if(isslr){
                                 analyses_temp[numanalyses++]=SealevelriseAnalysisEnum;
+                        }
+                        if(isesa){
+                                analyses_temp[numanalyses++]=EsaAnalysisEnum;
                         }
                         if(isgia){
@@ -2291 +2299 @@
 /*}}}*/
 #endif
+#ifdef _HAVE_ESA_
+void FemModel::EsaGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
+
+        IssmDouble  eartharea=0;
+        IssmDouble  eartharea_cpu=0;
+
+        int         ns,nsmax;
+        
+        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
+        ns = elements->Size();
+        
+        /*First, figure out the surface area of Earth: */ 
+        for(int i=0;i<ns;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());
+
+        /*Figure out max of ns: */
+        ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        /*Call the esa geodetic core: */
+        for(int i=0;i<nsmax;i++){
+                if(i<ns){
+                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                        element->EsaGeodetic(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea);
+                }
+                if(i%100==0){
+                        pUp->Assemble();
+                        pNorth->Assemble();
+                        pEast->Assemble();
+                }
+        }
+        
+        /*One last time: */
+        pUp->Assemble();
+        pNorth->Assemble();
+        pEast->Assemble();
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(latitude);
+        xDelete<IssmDouble>(longitude);
+        xDelete<IssmDouble>(radius);
+        xDelete<IssmDouble>(xx);
+        xDelete<IssmDouble>(yy);
+        xDelete<IssmDouble>(zz);
+}
+/*}}}*/
+#endif
 #ifdef _HAVE_SEALEVELRISE_
 void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/