Changeset 21395

Timestamp:

11/17/16 22:44:51 (8 years ago)

Author:

Eric.Larour

Message:

CHG: rotational feedback.

Location:

issm/branches/trunk-larour-NatClimateChange2016/src/c/classes

Files:

• FemModel.cpp (modified) (1 diff)
• FemModel.h (modified) (1 diff)

issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.cpp

@@ -2402 +2402 @@
 }
 /*}}}*/
+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, radi, 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];     // term associated with fluid number (3-order-of-magnitude smaller) is negelected  
+
+        /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
+         * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) 
+         * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) 
+         * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] 
+         */
+        for(int i=0;i<vertices->Size();i++){
+                int sid;
+                //Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
+                Vertex* vertex=xDynamicCast<Vertex*>(vertices->GetObjectByOffset(i));
+                sid=vertex->Sid();
+
+                lati=latitude[sid]/180*PI;      longi=longitude[sid]/180*PI; radi=radius[sid];
+
+                /*only first order terms are considered now: */ 
+                value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radi,2.0)*
+                                                (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); 
+        
+                pSgo_rot->SetValue(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);
+        
+}
+/*}}}*/
 void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
 

issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/FemModel.h

@@ -116 +116 @@
                 void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
                 void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
+                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
                 void SealevelriseGeodetic(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); 
                 IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);