Changeset 14807

Timestamp:

04/30/13 10:00:25 (12 years ago)

Author:

Eric.Larour

Message:

CHG: gia_core now return W and dW/dt
Added runtime controls for output_rates and cross_section_shape, which map into irate and iedge in naruse.

Location:

issm/trunk-jpl/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/classes/FemModel.cpp (modified) (2 diffs)
• c/classes/FemModel.h (modified) (1 diff)
• c/classes/GiaDeflectionCoreArgs.h (modified) (1 diff)
• c/classes/objects/Elements/Element.h (modified) (1 diff)
• c/classes/objects/Elements/Penta.cpp (modified) (2 diffs)
• c/classes/objects/Elements/Penta.h (modified) (1 diff)
• c/classes/objects/Elements/Tria.cpp (modified) (7 diffs)
• c/classes/objects/Elements/Tria.h (modified) (1 diff)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (modified) (4 diffs)
• c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (modified) (1 diff)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (2 diffs)
• c/solutions/gia_core.cpp (modified) (3 diffs)
• m/classes/gia.m (modified) (4 diffs)
• m/enum/EnumDefinitions.py (modified) (2 diffs)
• m/enum/GiaCrossSectionShapeEnum.m (added)
• m/enum/GiaOutputRatesEnum.m (added)
• m/enum/GiadWdtEnum.m (added)
• m/enum/MaximumNumberOfEnums.m (modified) (1 diff)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h ¶

@@ -477 +477 @@
         StressTensoryzEnum,
         StressTensorzzEnum,
+        GiaOutputRatesEnum,
+        GiaCrossSectionShapeEnum,
+        GiadWdtEnum,
         GiaWEnum,
         /*}}}*/

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

@@ -1544 +1544 @@
 #endif
 #ifdef _HAVE_GIA_
-void FemModel::Deflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){ /*{{{*/
+void FemModel::Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/
 
         int      i;
@@ -1554 +1554 @@
         for (i=0;i<elements->Size();i++){
                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                element->GiaDeflection(wg,x,y);
-        }
-}
-/*}}}*/
-#endif
+                element->GiaDeflection(wg,dwgdt, x,y);
+        }
+}
+/*}}}*/
+#endif

issm/trunk-jpl/src/c/classes/FemModel.h ¶

@@ -94 +94 @@
                 #endif
                 #ifdef _HAVE_GIA_
-                void Deflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y);
+                void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
                 #endif
                 void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax);

issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h ¶

@@ -25 +25 @@
         IssmDouble lithosphere_thickness;
 
+        /*gia solution parameters: */
+        int iedge;
+        int irate;
+
         /*ice properties: */
         IssmDouble rho_ice;

issm/trunk-jpl/src/c/classes/objects/Elements/Element.h ¶

@@ -101 +101 @@
 
                 #ifdef _HAVE_GIA_
-                virtual void   GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y)=0;
+                virtual void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y)=0;
                 #endif
 

issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp ¶

@@ -2171 +2171 @@
                                 name==QmuSurfaceEnum ||
                                 name==QmuTemperatureEnum ||
-                                name==QmuMeltingEnum
+                                name==QmuMeltingEnum ||
+                                name==GiaWEnum || 
+                                name==GiadWdtEnum
+
                                 ) {
                 return true;
@@ -3498 +3501 @@
 #ifdef _HAVE_GIA_
 /*FUNCTION Penta::GiaDeflection {{{*/
-void Penta::GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){
+void Penta::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y){
         _error_("GIA deflection not implemented yet!");
 }

issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h ¶

@@ -144 +144 @@
 
                 #ifdef _HAVE_GIA_
-                void   GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y);
+                void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
                 #endif
 

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

@@ -2060 +2060 @@
                                 name==QmuSurfaceEnum ||
                                 name==QmuTemperatureEnum ||
-                                name==QmuMeltingEnum
+                                name==QmuMeltingEnum || 
+                                name==GiaWEnum || 
+                                name==GiadWdtEnum
                 ){
                 return true;
@@ -3072 +3074 @@
 #ifdef _HAVE_GIA_
 /*FUNCTION Tria::GiaDeflection {{{*/
-void Tria::GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){
+void Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x, IssmDouble* y){
 
         int i;
@@ -3087 +3089 @@
         int         numtimes;
         Input* thickness_input=NULL;
+
+        /*gia solution parameters:*/
+        int output_rates=0;
+        int cross_section_shape=0;
 
         /*gia material parameters: */
@@ -3102 +3108 @@
         /*output: */
         IssmDouble  wi;
+        IssmDouble  dwidt;
 
         /*arguments to GiaDeflectionCorex: */
@@ -3108 +3115 @@
         /*how many dofs are we working with here? */
         this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
+
+        /*recover gia solution parameters: */
+        this->parameters->FindParam(&output_rates,GiaOutputRatesEnum);
+        this->parameters->FindParam(&cross_section_shape,GiaCrossSectionShapeEnum);
 
         /*what time is it? :*/
@@ -3157 +3168 @@
         arguments.rho_ice=rho_ice;
         arguments.idisk=this->id;
+        arguments.irate=output_rates;
+        arguments.iedge=cross_section_shape;
 
         for(i=0;i<gsize;i++){
@@ -3167 +3180 @@
 
                 /*for this Tria, compute contribution to rebound at vertex i: */
-                GiaDeflectionCorex(&wi,&arguments);
+                GiaDeflectionCorex(&wi,&dwidt,&arguments);
 
                 /*plug value into solution vector: */
                 wg->SetValue(i,wi,ADD_VAL);
+                dwgdt->SetValue(i,dwidt,ADD_VAL);
 
         }

issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h ¶

@@ -142 +142 @@
 
                 #ifdef _HAVE_GIA_
-                void   GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y);
+                void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
                 void   InputUpdateFromSolutionGia(IssmDouble* solution);
                 #endif

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp ¶

@@ -468 +468 @@
                 case StressTensoryzEnum : return "StressTensoryz";
                 case StressTensorzzEnum : return "StressTensorzz";
+                case GiaOutputRatesEnum : return "GiaOutputRates";
+                case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape";
+                case GiadWdtEnum : return "GiadWdt";
                 case GiaWEnum : return "GiaW";
                 case P0Enum : return "P0";

issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp ¶

@@ -68 +68 @@
 /*}}}*/
       
-void GiaDeflectionCorex( IssmDouble* pwi, GiaDeflectionCoreArgs* arguments){
+void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){
 
         /*output: */
         IssmDouble wi=0;
+        IssmDouble dwidt=0;
         
-        /*inputs: */
-        int iedge=1; //c iedge ......... = 1 square-edged, = 2 elliptical disc x-section  (see naruse.f)
+        /*inputs: {{{*/
+        /*constant: */
         int idisk=1; // disk #
-        int irate=0; // =0 fetch w solution (m) only; =1 dw/dt (mm/yr) only 
 
-        /*intermediary: */
+        /*coming from the model structure, runtime configurable:*/
+        /*gia solution parameters: */
+        int iedge=0; 
+        int irate=0;
 
-        /*inputs: {{{*/
+        /*gia inputs: */
         IssmDouble ri; //radial distance from center of disk to vertex  i
         IssmDouble re; //radius of disk
@@ -121 +124 @@
         rho_ice=arguments->rho_ice; 
         disk_id=arguments->idisk;
+        irate=arguments->irate;
+        iedge=arguments->iedge;
+
+        //printf("%g %g %g %i %g %g %g %g %g %g %g %g %i\n", ri,re,current_he,numtimes,currenttime,lithosphere_shear_modulus,lithosphere_density,mantle_shear_modulus,mantle_viscosity, mantle_density,lithosphere_thickness,rho_ice,disk_id);
+
         /*}}}*/
 
+        /*Modify inputs to match naruse code: */
+        //from our model, irate comes in with values in [1,2], which maps into [0,1] in naruse: 
+        irate=irate-1;
+
+        /*Prepare block inputs for fortran distme and what0 routines of the naruse code: {{{*/
         /*Now, let's set pset from the data that we got in input to GiaDeflectionCorex: */
         blockp_.pset[0]=lithosphere_thickness;
@@ -147 +160 @@
         /*radial distance of i-th element: */
         blockrad_.distrad=ri/1000.0; // in km
+        /*}}}*/
 
+        
         /*Call distme driver: */
         distme_(&idisk,&iedge); 
@@ -154 +169 @@
         what0_(&idisk,&iedge); 
 
-        /*this is the solution: */ 
-        wi = blocks_.aswokm;
+        /*output solution: */
+        if(irate==0){
+                wi = blocks_.aswokm;
+                *pwi=wi;
+                *pdwidt=0;
+        }
+        else if (irate==1){
+                dwidt = blocks_.aswokm;
+                *pdwidt=dwidt;
+                *pwi=0;
+        }
 
-        /*allocate output pointer: */
-        *pwi=wi;
-
-//      printf("wi: %g deflection: %g \n",wi,blocks_.aswokm); 
+        //printf("deflection: %g deflection rate %g\n",wi, dwidt);
 
 }

issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h ¶

@@ -10 +10 @@
 
 /* local prototypes: */
-void GiaDeflectionCorex( IssmDouble* pwi, GiaDeflectionCoreArgs* arguments);
+void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments);
 
 #endif  /* _GIADEFLECTIONCOREX_H */

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -105 +105 @@
         parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIssmbgradientsEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(GiaOutputRatesEnum));
 
         iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp ¶

@@ -478 +478 @@
               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
               else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+              else if (strcmp(name,"GiaOutputRates")==0) return GiaOutputRatesEnum;
+              else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
+              else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
               else if (strcmp(name,"GiaW")==0) return GiaWEnum;
               else if (strcmp(name,"P0")==0) return P0Enum;
@@ -504 +507 @@
               else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
               else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
-              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
-              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
-              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"MinVz")==0) return MinVzEnum;
+              if (strcmp(name,"MinVy")==0) return MinVyEnum;
+              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
+              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
+              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
               else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
               else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;

issm/trunk-jpl/src/c/solutions/gia_core.cpp ¶

@@ -15 +15 @@
         int i;
         Vector<IssmDouble>*  wg  = NULL;
+        Vector<IssmDouble>*  dwgdt  = NULL;
         IssmDouble*          x   = NULL;
         IssmDouble*          y   = NULL;
@@ -35 +36 @@
         gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
         wg = new Vector<IssmDouble>(gsize);
+        dwgdt = new Vector<IssmDouble>(gsize);
 
         /*first, recover x and y vectors from vertices: */
@@ -40 +42 @@
 
         /*call the main module: */
-        femmodel->Deflection(wg,x,y);
+        femmodel->Deflection(wg,dwgdt,x,y);
 
         /*assemble vector: */
         wg->Assemble();
+        dwgdt->Assemble();
 
         InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,wg); 
+        InputUpdateFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,dwgdt,GiadWdtEnum,VertexEnum); 
 
         if(save_results){
                 if(VerboseSolution()) _pprintLine_("   saving results");
                 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GiaWEnum);
+                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GiadWdtEnum);
         }
 

issm/trunk-jpl/src/m/classes/gia.m ¶

@@ -7 +7 @@
         properties (SetAccess=public) 
                 lithosphere_thickness         = NaN;
+                output_rates                  = 0;
+                cross_section_shape           = 0;
         end
         methods
@@ -18 +20 @@
                 end % }}}
                 function obj = setdefaultparameters(obj) % {{{
-
+                obj.output_rates=1;  %don't output  rates by default (see irate in GiaDeflectionCorex)
+                obj.cross_section_shape=1; %square as default (see iedge in GiaDeflectionCorex)
                 end % }}}
                 function md = checkconsistency(obj,md,solution,analyses) % {{{
@@ -24 +27 @@
                         if ~ismember(GiaAnalysisEnum(),analyses), return; end
                         md = checkfield(md,'gia.lithosphere_thickness','NaN',1,'size',[md.mesh.numberofvertices 1],'>',0);
+                        md = checkfield(md,'gia.output_rates','numel',[1],'values',[1,2]);
+                        md = checkfield(md,'gia.cross_section_shape','numel',[1],'values',[1,2]);
 
                 end % }}}
@@ -30 +35 @@
 
                         fielddisplay(obj,'lithosphere_thickness','lithosphere thickness[km]');
+                        fielddisplay(obj,'output_rates','1: fetch w solution (m) (default). 2: fetch dw/dt (mm/yr). See irate in GiaDeflectionCore');
+                        fielddisplay(obj,'cross_section_shape','1: square-edged (default). 2: elliptical.  See iedge in GiaDeflectionCore');
 
                 end % }}}
                 function marshall(obj,fid) % {{{
                         WriteData(fid,'data',obj.lithosphere_thickness,'format','DoubleMat','mattype',1,'enum',GiaLithosphereThicknessEnum());
+                        WriteData(fid,'object',obj,'fieldname','output_rates','format','Integer');
+                        WriteData(fid,'object',obj,'fieldname','cross_section_shape','format','Integer');
                 end % }}}
         end

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

@@ -4519 +4519 @@
         return StringToEnum('StressTensorzz')[0]
 
+def GiaOutputRatesEnum():
+        """
+        GIAOUTPUTRATESENUM - Enum of GiaOutputRates
+
+           Usage:
+              macro=GiaOutputRatesEnum()
+        """
+
+        return StringToEnum('GiaOutputRates')[0]
+
+def GiaCrossSectionShapeEnum():
+        """
+        GIACROSSSECTIONSHAPEENUM - Enum of GiaCrossSectionShape
+
+           Usage:
+              macro=GiaCrossSectionShapeEnum()
+        """
+
+        return StringToEnum('GiaCrossSectionShape')[0]
+
+def GiadWdtEnum():
+        """
+        GIADWDTENUM - Enum of GiadWdt
+
+           Usage:
+              macro=GiadWdtEnum()
+        """
+
+        return StringToEnum('GiadWdt')[0]
+
 def GiaWEnum():
         """
@@ -5427 +5457 @@
         """
 
-        return 541
-
+        return 544
+

issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m ¶

@@ -9 +9 @@
 %      macro=MaximumNumberOfEnums()
 
-macro=541;
+macro=544;