Changeset 18551

Timestamp:

09/30/14 11:03:47 (10 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on balancethickness2 (added new field), and some minor changes to googlemaps to make it work for other regions than GIS or AIS

Location:

issm/trunk-jpl/src

Files:

• c/analyses/Balancethickness2Analysis.cpp (modified) (7 diffs)
• c/cores/balancethickness2_core.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• m/classes/balancethickness.m (modified) (3 diffs)
• m/enum/BalancethicknessCmuEnum.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/plot/googlemaps.m (modified) (1 diff)
• m/plot/plot_googlemaps.m (modified) (1 diff)

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

@@ -34 +34 @@
         iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
         iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
+        iomodel->FetchDataToInput(elements,BalancethicknessCmuEnum);
 }/*}}}*/
 void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -144 +145 @@
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
                                         (ms-mb-dhdt)*basis[i]
-                                        -D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
+                                        //-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
                                         );
         }
@@ -155 +156 @@
 }/*}}}*/
 void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-           _error_("not implemented yet");
+                element->GetSolutionFromInputsOneDof(solution,SurfaceEnum);
 }/*}}}*/
 void Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -161 +162 @@
 }/*}}}*/
 void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
+
+        /*Intermediaries*/
+        IssmDouble  ds[2],s,b,D;
+        IssmDouble* xyz_list = NULL;
+
+        //element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
+        element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
+
+        /*Fetch number of vertices and allocate velocity vectors*/
+        int numvertices = element->GetNumberOfVertices();
+        IssmDouble* vel_list = xNew<IssmDouble>(numvertices);
+        IssmDouble* vx_list  = xNew<IssmDouble>(numvertices);
+        IssmDouble* vy_list  = xNew<IssmDouble>(numvertices);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);
+        Input* H_input   = element->GetInput(ThicknessEnum);                            _assert_(H_input);
+        Input* s_input   = element->GetInput(SurfaceEnum);                              _assert_(s_input);
+        Input* b_input   = element->GetInput(BaseEnum);                                 _assert_(b_input);
+
+        /*Calculate velocities*/
+        Gauss* gauss=element->NewGauss();
+        for(int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+
+                if(D_input){
+                        D_input->GetInputValue(&D,gauss);
+                }
+                else{
+                        D = 0.;
+                }
+                b_input->GetInputValue(&b,gauss);
+                s_input->GetInputValue(&s,gauss);
+                s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+
+                vx_list[iv] = -1./(s-b)*D*ds[0];
+                vy_list[iv] = -1./(s-b)*D*ds[1];
+                vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2));
+        }
+
+        /*Add vx and vy as inputs to the tria element: */
+        element->AddInput(VxEnum,vx_list,P1Enum);
+        element->AddInput(VyEnum,vy_list,P1Enum);
+        element->AddInput(VelEnum,vel_list,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(vy_list);
+        xDelete<IssmDouble>(vx_list);
+        xDelete<IssmDouble>(vel_list);
+        xDelete<IssmDouble>(xyz_list);
 }/*}}}*/
 void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -172 +223 @@
 
         /*Intermediaries */
-        IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k;
+        IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b;
         const int        n = 3;
         const IssmDouble Hstar = 500.;
@@ -185 +236 @@
         /*retrieve what we need: */
         element->GetVerticesCoordinates(&xyz_list);
-        Input* thickness_input  = element->GetInput(ThicknessEnum);          _assert_(thickness_input);
+        Input* bed_input        = element->GetInput(BaseEnum);               _assert_(bed_input);
         Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
         Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
         Input* k_input          = element->GetInput(FrictionCoefficientEnum);_assert_(k_input);
+        Input* Cmu_input        = element->GetInput(BalancethicknessCmuEnum);_assert_(Cmu_input);
         IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
@@ -198 +250 @@
                 B_input->GetInputValue(&B,gauss);
                 k_input->GetInputValue(&k,gauss);
-                thickness_input->GetInputValue(&h,gauss);
+                //thickness_input->GetInputValue(&h,gauss);
                 surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+                surface_input->GetInputValue(&s,gauss);
+                bed_input->GetInputValue(&b,gauss);
+                h = s-b;
 
                 mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
                 Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
-                Cmu   = mu0/Hstar/(k*k);
-
-                D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h - Cmu);
+                //Cmu   = pow(mu0,n)/Hstar/(k*k);
+                Cmu_input->GetInputValue(&Cmu,gauss);
+
+                D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h + Cmu);
         }
 

issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp

@@ -22 +22 @@
 
         if(VerboseSolution()) _printf0_("call computational core:\n");
-        solutionsequence_linear(femmodel);
+        //solutionsequence_linear(femmodel);
+        solutionsequence_nonlinear(femmodel,false);
 
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                int outputs[1] = {ThicknessEnum};
-                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
+                const int numoutputs = 5;
+                //int outputs[numoutputs] = {ThicknessEnum};
+                int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum,BalancethicknessDiffusionCoefficientEnum};
+                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
         }
 

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

@@ -305 +305 @@
         Balancethickness2MisfitEnum,
         BalancethicknessDiffusionCoefficientEnum,
+        BalancethicknessCmuEnum,
         /*}}}*/
         /*Surfaceforcings{{{*/

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

@@ -313 +313 @@
                 case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
                 case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient";
+                case BalancethicknessCmuEnum : return "BalancethicknessCmu";
                 case SurfaceforcingsEnum : return "Surfaceforcings";
                 case SMBEnum : return "SMB";

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

@@ -319 +319 @@
               else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum;
               else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
+              else if (strcmp(name,"BalancethicknessCmu")==0) return BalancethicknessCmuEnum;
               else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum;
               else if (strcmp(name,"SMB")==0) return SMBEnum;
@@ -382 +383 @@
               else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
               else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
-              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
+              if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
               else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
               else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
@@ -505 +506 @@
               else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
               else if (strcmp(name,"Friction")==0) return FrictionEnum;
-              else if (strcmp(name,"Internal")==0) return InternalEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
+              if (strcmp(name,"Internal")==0) return InternalEnum;
+              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
               else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
               else if (strcmp(name,"Misfit")==0) return MisfitEnum;
@@ -628 +629 @@
               else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
               else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
-              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
+              if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+              else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
               else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
               else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
@@ -751 +752 @@
               else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
               else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
-              else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
+              if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
+              else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
               else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
               else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;

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

@@ -9 +9 @@
                 thickening_rate   = NaN;
                 stabilization     = 0;
+
+                Cmu               = NaN;
         end
         methods
-        function createxml(obj,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- balancethickness -->');
-            
-                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcthickness" type="',     class(obj.spcthickness),'" default="',        convert2str(obj.spcthickness),'">',       '     <section name="balancethickness" />','     <help> thickness constraints (NaN means no constraint) [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="thickening_rate" type="',      class(obj.thickening_rate),'" default="',            convert2str(obj.thickening_rate),'">',     '     <section name="balancethickness" />','     <help> ice thickening rate used in the mass conservation (dh/dt) [m/yr] </help>','</parameter>');
-            
-            % balancethickness drop-down (1,2, or 3)
-            fprintf(fid,'%s\n%s\n%s\n',        '<parameter key ="stabilization" type="alternative" optional="false">','     <section name="balancethickness" />','     <help> 0: None, 1: SU, 2: SSAs artificial diffusivity, 3:DG </help>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n','       <option value="2" type="string" default="false"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="3" type="string" default="false"> </option>','</parameter>');
-
-        end % }}}
                 function obj = balancethickness(varargin) % {{{
                         switch nargin
@@ -46 +34 @@
                         md = checkfield(md,'fieldname','balancethickness.thickening_rate','size',[md.mesh.numberofvertices 1],'NaN',1);
                         md = checkfield(md,'fieldname','balancethickness.stabilization','size',[1 1],'values',[0 1 2 3]);
+
+                        md = checkfield(md,'fieldname','balancethickness.Cmu','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0);
                 end % }}}
                 function disp(obj) % {{{
@@ -62 +52 @@
                         WriteData(fid,'object',obj,'fieldname','thickening_rate','format','DoubleMat','mattype',1,'scale',1./yts);
                         WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
+
+                        WriteData(fid,'object',obj,'fieldname','Cmu','format','DoubleMat','mattype',1);
                 end % }}}
         end

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -305 +305 @@
 def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0]
 def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0]
+def BalancethicknessCmuEnum(): return StringToEnum("BalancethicknessCmu")[0]
 def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0]
 def SMBEnum(): return StringToEnum("SMB")[0]

issm/trunk-jpl/src/m/plot/googlemaps.m

@@ -54 +54 @@
 elseif strcmpi(md.mesh.hemisphere,'s'),
         EPSGlocal = 'EPSG:3031'; % UPS Antarctica http://www.spatialreference.org/ref/epsg/3031/
+elseif strcmpi(md.mesh.hemisphere,'UTM18N'),
+        EPSGlocal = 'EPSG:32218'; % UTM 18 N Patagonia
 else
         error('field hemisphere should either be ''n'' or ''s''');

issm/trunk-jpl/src/m/plot/plot_googlemaps.m

@@ -36 +36 @@
                         -1,0,71);
         else
+                latlist = md.mesh.lat; %That might work?
+                lonlist = md.mesh.long;
                 error('field hemisphere should either be ''n'' or ''s''');
         end