Changeset 26232

Timestamp:

05/04/21 09:11:16 (4 years ago)

Author:

caronlam

Message:

Adding support for rotationnal feedback in bed Up, East, North

File:

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

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

@@ -6133 +6133 @@
         IssmDouble* tide_love_h  = NULL;
         IssmDouble* tide_love_k  = NULL;
-        IssmDouble* load_love_k  = NULL;
-        IssmDouble  tide_love_k2secular;
+        IssmDouble* tide_love_l  = NULL;
         IssmDouble  moi_e, moi_p, omega;
-        IssmDouble Grotm1[3];
-        IssmDouble Grotm2[3];
-        IssmDouble Grotm3[3];
-        IssmDouble pre;
+        IssmDouble  Grotm1[3],GUrotm1[3],GNrotm1[3],GErotm1[3];
+        IssmDouble  Grotm2[3],GUrotm2[3],GNrotm2[3],GErotm2[3];
+        IssmDouble  Grotm3[3],GUrotm3[3],GNrotm3[3],GErotm3[3];
+        IssmDouble  Y21cos     , Y21sin     , Y20;
+        IssmDouble dY21cos_dlat,dY21sin_dlat,dY20_dlat;
+        IssmDouble dY21cos_dlon,dY21sin_dlon;
+        IssmDouble LoveRotRSL,LoveRotU,LoveRothoriz;
         /*}}}*/
         /*Recover parameters:{{{ */
@@ -6153 +6155 @@
 
         if(computerotation){
-                parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
                 parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
                 parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
-                parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
+                parameters->FindParam(&tide_love_l,NULL,NULL,TidalLoveLEnum);
                 parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
                 parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
@@ -6270 +6271 @@
                 g=4.0/3.0*M_PI*rho_earth*NewtonG*planetradius;
 
+                //Amplitude of the rotational feedback
+                LoveRotRSL=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0);
+                LoveRotU=(tide_love_h[2]/g)*pow(omega*planetradius,2.0);
+                LoveRothoriz=(tide_love_l[2]/g)*pow(omega*planetradius,2.0);
+
                 for(int i=0;i<3;i++){
                         lati=latitude[i];
                         longi=longitude[i];
 
-                        pre=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0);
-                        Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi);
-                        Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi);
-                        Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati));
+                        //Spherical harmonic functions of degree 2
+                        Y21cos= -0.5*sin(2.*lati)*cos(longi);
+                        Y21sin= -0.5*sin(2.*lati)*sin(longi);
+                        Y20   = -(1.0/6.0 - 0.5*cos(2.0*lati));
+
+                        Grotm1[i]=  LoveRotRSL*Y21cos;
+                        Grotm2[i]=  LoveRotRSL*Y21sin;
+                        Grotm3[i]=  LoveRotRSL*Y20;
+
+                        if (computeelastic){
+                                GUrotm1[i]=  LoveRotU*Y21cos;
+                                GUrotm2[i]=  LoveRotU*Y21sin;
+                                GUrotm3[i]=  LoveRotU*Y20;
+                                if (horiz){
+                                //bed_N = Love_l * d(Y)/dlat ;
+                                dY21cos_dlat=-cos(2.*lati)*cos(longi);
+                                dY21sin_dlat=-cos(2.*lati)*sin(longi);
+                                dY20_dlat=-sin(2.0*lati);
+                                GNrotm1[i]=  LoveRothoriz*dY21cos_dlat;
+                                GNrotm2[i]=  LoveRothoriz*dY21sin_dlat;
+                                GNrotm3[i]=  LoveRothoriz*dY20_dlat;
+
+                                //bed_E = Love_l * 1/cos(lat) * d(Y)/dlon ;
+                                dY21cos_dlon=-Y21sin;
+                                dY21sin_dlon=Y21cos;
+                                //dY20_dlon=0.;
+                                GNrotm1[i]=  LoveRothoriz*dY21cos_dlon/cos(lati);
+                                GNrotm2[i]=  LoveRothoriz*dY21sin_dlon/cos(lati);
+                                GNrotm3[i]=  0.0;
+                                }
+                        }
                 }
                 this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
                 this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum);
                 this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum);
+                if (computeelastic){
+                        this->AddInput(SealevelGUrotm1Enum,&GUrotm1[0],P1Enum);
+                        this->AddInput(SealevelGUrotm2Enum,&GUrotm2[0],P1Enum);
+                        this->AddInput(SealevelGUrotm3Enum,&GUrotm3[0],P1Enum);
+                        if(horiz){
+                                this->AddInput(SealevelGNrotm1Enum,&GNrotm1[0],P1Enum);
+                                this->AddInput(SealevelGNrotm2Enum,&GNrotm2[0],P1Enum);
+                                this->AddInput(SealevelGNrotm3Enum,&GNrotm3[0],P1Enum);
+                                this->AddInput(SealevelGErotm1Enum,&GErotm1[0],P1Enum);
+                                this->AddInput(SealevelGErotm2Enum,&GErotm2[0],P1Enum);
+                                this->AddInput(SealevelGErotm3Enum,&GErotm3[0],P1Enum);
+                        }
+                }
         }
         /*}}}*/
@@ -6995 +7041 @@
         IssmDouble Grotm2[3];
         IssmDouble Grotm3[3];
+        IssmDouble GUrotm1[3];
+        IssmDouble GUrotm2[3];
+        IssmDouble GUrotm3[3];
+        IssmDouble GNrotm1[3];
+        IssmDouble GNrotm2[3];
+        IssmDouble GNrotm3[3];
+        IssmDouble GErotm1[3];
+        IssmDouble GErotm2[3];
+        IssmDouble GErotm3[3];
         bool rotation= false;
         bool elastic=false;
@@ -7053 +7108 @@
                         }
                 }
+
+                if (elastic){
+                        Element::GetInputListOnVertices(&GUrotm1[0],SealevelGUrotm1Enum);
+                        Element::GetInputListOnVertices(&GUrotm2[0],SealevelGUrotm2Enum);
+                        Element::GetInputListOnVertices(&GUrotm3[0],SealevelGUrotm3Enum);
+                
+                        for(int i=0;i<NUMVERTICES;i++){
+                                if(slgeom->lids[this->vertices[i]->lid]==this->lid){
+                                        SealevelU[i]+=GUrotm1[i]*rotationvector[0]+GUrotm2[i]*rotationvector[1]+GUrotm3[i]*rotationvector[2];
+                                }
+                        }
+                        if (horiz){
+                                Element::GetInputListOnVertices(&GNrotm1[0],SealevelGNrotm1Enum);
+                                Element::GetInputListOnVertices(&GNrotm2[0],SealevelGNrotm2Enum);
+                                Element::GetInputListOnVertices(&GNrotm3[0],SealevelGNrotm3Enum);
+                                Element::GetInputListOnVertices(&GErotm1[0],SealevelGErotm1Enum);
+                                Element::GetInputListOnVertices(&GErotm2[0],SealevelGErotm2Enum);
+                                Element::GetInputListOnVertices(&GErotm3[0],SealevelGErotm3Enum);
+                
+                                for(int i=0;i<NUMVERTICES;i++){
+                                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
+                                                SealevelN[i]+=GNrotm1[i]*rotationvector[0]+GNrotm2[i]*rotationvector[1]+GNrotm3[i]*rotationvector[2];
+                                                SealevelE[i]+=GErotm1[i]*rotationvector[0]+GErotm2[i]*rotationvector[1]+GErotm3[i]*rotationvector[2];
+                                        }
+                                }
+                        }
+                }
         }
 
@@ -7062 +7144 @@
         this->AddInput(BedGRDEnum,SealevelU,P1Enum);
         if(horiz){
+                this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum);
                 this->AddInput(BedEastGRDEnum,SealevelE,P1Enum);
-                this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum);
         }