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

Timestamp:

11/04/16 17:42:57 (8 years ago)

Author:

adhikari

Message:

CHG: rotational feedback works fine now

File:

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

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

@@ -2459 +2459 @@
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-        
+
         /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
         for(int i=0;i<ns;i++){
@@ -2465 +2465 @@
                 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());
@@ -2475 +2476 @@
         for(int i=0;i<nsmax;i++){
                 if(i<ns){
-
                         if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
-
                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,eartharea);
+                        element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
                 }
                 if(i%100==0)pSgo->Assemble();
         }
         if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-
+        
         /*Free ressources:*/
         xDelete<IssmDouble>(Sg_old);
+}
+/*}}}*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+
+        /*serialized vectors:*/
+        IssmDouble* Sg_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 
+        IssmDouble      m1, m2, m3; 
+        IssmDouble      lati, longi, value; 
+
+        /*Serialize vectors from previous iteration:*/
+        Sg_old=pSg_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}; 
+        IssmDouble moi_list_cpu[3]={0,0,0}; 
+        for(int i=0;i<elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
+                moi_list_cpu[0] += moi_list[0]; 
+                moi_list_cpu[1] += moi_list[1]; 
+                moi_list_cpu[2] += moi_list[2]; 
+        }
+        ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        //      
+        ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        //      
+        ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        
+        /*pull out some useful parameters: */
+        parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum);
+        parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
+        parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum);
+        parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum);
+        parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum);
+        parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
+
+        /*compute perturbation terms for angular velocity vector: */
+        m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; 
+        m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; 
+        m3 = -(1+load_love_k2)/moi_p * moi_list[2]; 
+
+        for(int i=0;i<vertices->Size();i++){
+                //Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
+                Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
+                
+                lati=latitude[i]/180*PI;        longi=longitude[i]/180*PI;
+
+                /*only first order terms are considered now: */ 
+                value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radius[i],2.0)*
+                                                (-m3/6.0 - m3*cos(2.0*lati)/2.0 - sin(lati)*cos(lati)*(m1*cos(longi)+m2*sin(longi))); 
+        
+                pSgo_rot->SetValue(vertex->Sid(),value,INS_VAL); //INS_VAL ensures that you don't add several times
+        }
+
+        /*Assemble mesh velocity*/
+        pSgo_rot->Assemble();
+        
+        /*Free ressources:*/
+        xDelete<IssmDouble>(Sg_old);
+        
 }
 /*}}}*/