Changeset 21331

Timestamp:

11/01/16 11:48:51 (8 years ago)

Author:

adhikari

Message:

CHG: rotational feedback added (have to define perturbation terms for intertia tensor though)

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (4 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (2 diffs)
• m/classes/slr.js (modified) (7 diffs)
• m/classes/slr.m (modified) (7 diffs)
• m/classes/slr.py (modified) (6 diffs)

issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

@@ -70 +70 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
 

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

@@ -4048 +4048 @@
         bool computerigid = true;
         bool computeelastic= true;
+        bool computerotation= true;
 
         /*early return if we are not on the ocean:*/
@@ -4055 +4056 @@
         this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
         this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
-        
+        this->parameters->FindParam(&computerotation,SealevelriseRotationEnum);
+
+        /*recover some parameters for rotational feedback: */
+        IssmDouble tide_love_h, tide_love_k, omega;
+        if(computerotation){
+                this->parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum);
+                this->parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum);
+                this->parameters->FindParam(&omega,SealevelriseAngularVelocityEnum);
+        }
+
         /*early return if rigid or elastic not requested:*/
         if(!computerigid && !computeelastic) return;
@@ -4124 +4134 @@
         IssmDouble delPhi,delLambda;
 
+         /* if rotation, compute m1 m2 m3 from loading function */
+         IssmDouble m1, m2, m3;
+         if(computerotation){
+                 m1=0; 
+                 m2=0; 
+                 m3=0; 
+         }
+
         for(int i=0;i<gsize;i++){
 
@@ -4145 +4163 @@
                         G_elastic[i] += G_elastic_precomputed[index];
                         values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
+                }
+
+                /*Rotational feedback: */
+                if(computerotation){
+                        values[i]+=(1+tide_love_k-tide_love_h/9.81)*0.5*pow((omega*radius[i]),2)*(-m3*(2+m3)*(1+3*cos(2*lati))
+                                                + 2*pow(m2,2)*(pow(cos(lati),2) + 0.5*(-1 + 3*cos(2*longi)) * pow(sin(lati),2)) 
+                                                + 2*pow(m1,2)*(pow(cos(lati),2) - 0.5*(1+3*cos(2*longi)) * pow(sin(lati),2))
+                                                - m1*(1+m3)*cos(longi)*sin(2*lati)
+                       - 2*m1*m2*cos(longi)*pow(sin(lati),2)*sin(longi)
+                       - m2*(1+m3)*sin(2*lati)*sin(longi));
                 }
         }

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -773 +773 @@
         SealevelriseRigidEnum,
         SealevelriseElasticEnum,
+        SealevelriseRotationEnum,
+        SealevelriseTidalLoveHEnum,
+   SealevelriseTidalLoveKEnum,
+        SealevelriseAngularVelocityEnum,
         SealevelriseOceanAreaScalingEnum,
         SealevelriseGElasticEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -754 +754 @@
                 case SealevelriseRigidEnum : return "SealevelriseRigid";
                 case SealevelriseElasticEnum : return "SealevelriseElastic";
+                case SealevelriseRotationEnum : return "SealevelriseRotation";
+                case SealevelriseTidalLoveHEnum : return "SealevelriseTidalLoveH";
+                case SealevelriseTidalLoveKEnum : return "SealevelriseTidalLoveK";
+                case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity";
                 case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling";
                 case SealevelriseGElasticEnum : return "SealevelriseGElastic";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -772 +772 @@
               else if (strcmp(name,"SealevelriseRigid")==0) return SealevelriseRigidEnum;
               else if (strcmp(name,"SealevelriseElastic")==0) return SealevelriseElasticEnum;
+              else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum;
+              else if (strcmp(name,"SealevelriseTidalLoveH")==0) return SealevelriseTidalLoveHEnum;
+              else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum;
+              else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum;
               else if (strcmp(name,"SealevelriseOceanAreaScaling")==0) return SealevelriseOceanAreaScalingEnum;
               else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
@@ -871 +875 @@
               else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
               else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
-              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
+         else stage=8;
+   }
+   if(stage==8){
+              if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
               else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
               else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
               else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
-         else stage=8;
-   }
-   if(stage==8){
-              if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
+              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
               else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
               else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;

issm/trunk-jpl/src/m/classes/slr.js

@@ -18 +18 @@
                 this.rigid=1;
                 this.elastic=1;
-                this.rotation=1;
+                this.rotation=0;
                 this.ocean_area_scaling=0;
                 
@@ -24 +24 @@
                 this.tide_love_h=0.6149; //degree 2
                 this.tide_love_k=0.3055; //degree 2
+        
+                // mean rotational velocity of earth:
+                this.angular_velocity=7.2921*10^-5; // [s^-1] 
 
                 //numerical discretization accuracy
@@ -47 +50 @@
                         md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
+                        md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','slr.reltol','size',[1, 1]);
                         md = checkfield(md,'fieldname','slr.abstol','size',[1, 1]);
@@ -79 +83 @@
                 fielddisplay(this,'tide_love_h','tidal love number (degree 2)');
                 fielddisplay(this,'tide_love_k','tidal love number (degree 2)');
+                fielddisplay(this,'angular_velocity','mean rotational velocity of earth [per second]'); 
                 fielddisplay(this,'rigid','rigid earth graviational potential perturbation');
                 fielddisplay(this,'elastic','elastic earth graviational potential perturbation');
@@ -99 +104 @@
                         WriteData(fid,prefix,'object',this,'fieldname','tide_love_h','format','Double');
                         WriteData(fid,prefix,'object',this,'fieldname','tide_love_k','format','Double');
+                        WriteData(fid,prefix,'object',this,'fieldname','angular_velocity','format','Double');
                         WriteData(fid,prefix,'object',this,'fieldname','rigid','format','Boolean');
                         WriteData(fid,prefix,'object',this,'fieldname','elastic','format','Boolean');
@@ -128 +134 @@
                         this.tide_love_h=NullFix(this.tide_love_h,NaN);
                         this.tide_love_k=NullFix(this.tide_love_k,NaN);
+                        this.angular_velocity=NullFix(this.angular_velocity,NaN);
                         this.rigid=NullFix(this.rigid,NaN);
                         this.elastic=NullFix(this.elastic,NaN);
@@ -146 +153 @@
         this.tide_love_h    = 0; 
         this.tide_love_k    = 0; 
+        this.angular_velocity = 0; 
         this.rigid          = 0;
         this.elastic        = 0;

issm/trunk-jpl/src/m/classes/slr.m

@@ -16 +16 @@
                 tide_love_k    = 0; %ideam
                 tide_love_h    = 0; %ideam
+                angular_velocity = 0;
                 rigid          = 0;
                 elastic        = 0;
@@ -45 +46 @@
                 self.rigid=1;
                 self.elastic=1;
-                self.rotation=1;
+                self.rotation=0;
                 self.ocean_area_scaling=0;
 
@@ -51 +52 @@
                 self.tide_love_h=0.6149; %degree 2
                 self.tide_love_k=0.3055; % degree 2
+                
+                % mean rotational velocity of earth 
+                self.angular_velocity=7.2921*10^-5; % [s^-1] 
 
                 %numerical discretization accuracy
@@ -72 +76 @@
                         md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
+                        md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
                         md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
@@ -108 +113 @@
                         fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)');
                         fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)');
+                        fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'); 
                         fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
                         fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
@@ -128 +134 @@
                         WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
+                        WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean');
@@ -157 +164 @@
                         writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
                         writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
+                        writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
                         writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
                         writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);

issm/trunk-jpl/src/m/classes/slr.py

@@ -25 +25 @@
                 self.tide_love_h       = 0
                 self.tide_love_k       = 0
+                self.angular_velocity  = 0;
                 self.rigid             = 0
                 self.elastic           = 0
@@ -47 +48 @@
                         string="%s\n%s"%(string,fielddisplay(self,'tide_love_k','tidal load Love number (degree 2)'))
                         string="%s\n%s"%(string,fielddisplay(self,'tide_love_h','tidal load Love number (degree 2)'))
+                        string="%s\n%s"%(string,fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]')); 
                         string="%s\n%s"%(string,fielddisplay(self,'rigid','rigid earth graviational potential perturbation'))
                         string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation'))
@@ -69 +71 @@
                 self.rigid=1
                 self.elastic=1
-                self.rotation=1
+                self.rotation=0
                 self.ocean_area_scaling=0
 
@@ -75 +77 @@
                 self.tide_love_h=0.6149; #degree 2
                 self.tide_love_k=0.3055; #degree 2
+
+                #mean rotational velocity of earth 
+                self.angular_velocity=7.2921*10^-5; # [s^-1] 
 
                 #numerical discretization accuracy
@@ -102 +107 @@
                 md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1)
                 md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1)
+                md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1)
                 md = checkfield(md,'fieldname','slr.reltol','size',[1,1])
                 md = checkfield(md,'fieldname','slr.abstol','size',[1,1])
@@ -127 +133 @@
                 WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
                 WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
+                WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
                 WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean')
                 WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean')