Changeset 7848

Timestamp:

04/12/11 11:03:35 (14 years ago)

Author:

Mathieu Morlighem

Message:

Added rhelogy law for B temperature dependence
NoneEnum -> B does not depend on T
PatersonEnm -> use Paterson's table (default)
ArrheniusEnum -> use Arrhenius law from EISMINT tests

Location:

issm/trunk/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/EnumDefinitions/EnumToString.cpp (modified) (1 diff)
• c/EnumDefinitions/StringToEnum.cpp (modified) (1 diff)
• c/Makefile.am (modified) (2 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• c/objects/Elements/Penta.cpp (modified) (2 diffs)
• c/objects/IoModel.cpp (modified) (3 diffs)
• c/objects/IoModel.h (modified) (1 diff)
• c/shared/Elements/Arrhenius.cpp (added)
• c/shared/Elements/elements.h (modified) (1 diff)
• m/classes/model.m (modified) (2 diffs)
• m/enum/ArrheniusEnum.m (added)
• m/enum/EnumToString.m (modified) (1 diff)
• m/enum/PatersonEnum.m (added)
• m/enum/RheologyLawEnum.m (added)
• m/enum/StringToEnum.m (modified) (1 diff)
• m/model/marshall.m (modified) (1 diff)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -431 +431 @@
         OptionDoubleEnum,
         OptionLogicalEnum,
-        OptionStructEnum
-          /*}}}*/
+        OptionStructEnum,
+        /*}}}*/
+        /*Rheology law (move too Material) {{{1*/
+        RheologyLawEnum,
+        PatersonEnum,
+        ArrheniusEnum
+        /*}}}*/
 };
 

issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

@@ -381 +381 @@
                 case OptionLogicalEnum : return "OptionLogical";
                 case OptionStructEnum : return "OptionStruct";
+                case RheologyLawEnum : return "RheologyLaw";
+                case PatersonEnum : return "Paterson";
+                case ArrheniusEnum : return "Arrhenius";
                 default : return "unknown";
 

issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

@@ -379 +379 @@
         else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum;
         else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
+        else if (strcmp(name,"RheologyLaw")==0) return RheologyLawEnum;
+        else if (strcmp(name,"Paterson")==0) return PatersonEnum;
+        else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
         else _error_("Enum %s not found",name);
 

issm/trunk/src/c/Makefile.am

@@ -327 +327 @@
                                         ./shared/Elements/elements.h\
                                         ./shared/Elements/Paterson.cpp\
+                                        ./shared/Elements/Arrhenius.cpp\
                                         ./shared/Elements/GetVerticesCoordinates.cpp\
                                         ./shared/Elements/GetLocalDofList.cpp\
@@ -962 +963 @@
                                         ./shared/Elements/elements.h\
                                         ./shared/Elements/Paterson.cpp\
+                                        ./shared/Elements/Arrhenius.cpp\
                                         ./shared/Elements/GetVerticesCoordinates.cpp\
                                         ./shared/Elements/GetLocalDofList.cpp\

issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -36 +36 @@
         parameters = new Parameters(ParametersEnum);
 
-        parameters->AddObject(new    IntParam(SolutionTypeEnum,solution_type));
-        parameters->AddObject(new    IntParam(AnalysisTypeEnum,analysis_type));
-        parameters->AddObject(new    IntParam(AnalysisCounterEnum,analysis_counter));
-        if (iomodel->dim==2) parameters->AddObject(new IntParam(DimEnum,2));
-        else parameters->AddObject(new IntParam(DimEnum,3));
+        parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type));
+        parameters->AddObject(new IntParam(AnalysisTypeEnum,analysis_type));
+        parameters->AddObject(new IntParam(AnalysisCounterEnum,analysis_counter));
+
+        parameters->AddObject(new IntParam(DimEnum,iomodel->dim));
         parameters->AddObject(new BoolParam(IsHutterEnum,iomodel->ishutter));
         parameters->AddObject(new BoolParam(IsMacAyealPattynEnum,iomodel->ismacayealpattyn));
@@ -79 +79 @@
         parameters->AddObject(new BoolParam(IoGatherEnum,iomodel->io_gather));
         parameters->AddObject(new IntParam(GroundingLineMigrationEnum,iomodel->gl_migration));
+        parameters->AddObject(new IntParam(RheologyLawEnum,iomodel->rheology_law));
 
         /*Deal with more complex parameters*/

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -5095 +5095 @@
 
         bool   converged;
-        int    i;
+        int    i,rheology_law;
+        double xyz_list[NUMVERTICES][3];
         double values[numdof];
         double B[numdof];
-        double B_average;
+        double B_average,s_average;
         int*   doflist=NULL;
 
@@ -5114 +5115 @@
         }
 
+        /*Get all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+
         this->inputs->GetParameterValue(&converged,ConvergedEnum);
         if(converged){
                 this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,values));
 
-                /*Update Rheology only if convreged (we must make sure that the temperature is below melting point
+                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
                  * otherwise the rheology could be negative*/
-                B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
-                for(i=0;i<numdof;i++) B[i]=B_average;
-                this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
+                this->parameters->FindParam(&rheology_law,RheologyLawEnum);
+                switch(rheology_law){
+                        case NoneEnum:
+                                /*Do nothing: B is not temperature dependent*/
+                                break;
+                        case PatersonEnum:
+                                B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
+                                for(i=0;i<numdof;i++) B[i]=B_average;
+                                this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
+                                break;
+                        case ArrheniusEnum:
+                                surface_input->GetParameterAverage(&s_average);
+                                B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
+                                                        s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
+                                                        matice->GetN());
+                                for(i=0;i<numdof;i++) B[i]=B_average;
+                                _error_("Not supported yet (See Rutt2009 for example)");
+                                break;
+                        default:
+                                _error_("Rheology law %s not supported yet",EnumToString(rheology_law));
+
+                }
         }
         else{

issm/trunk/src/c/objects/IoModel.cpp

@@ -214 +214 @@
         IoModelFetchData(&this->melting_rate_correction_apply,iomodel_handle,"melting_rate_correction_apply");
         IoModelFetchData(&this->gl_melting_rate,iomodel_handle,"gl_melting_rate");
+        IoModelFetchData(&this->rheology_law,iomodel_handle,"rheology_law");
         
         /*qmu: */
@@ -225 +226 @@
                 IoModelFetchData(&this->qmu_save_femmodel,iomodel_handle,"qmu_save_femmodel"); 
         }
-
 
         /*i/o: */
@@ -321 +321 @@
         this->rheology_n=NULL;
         this->rheology_B=NULL;
+        this->rheology_law=0;
 
         /*!solution parameters: */

issm/trunk/src/c/objects/IoModel.h

@@ -113 +113 @@
                 double* rheology_B;
                 double* rheology_n;
+                int     rheology_law;
 
                 /*numerical parameters: */

issm/trunk/src/c/shared/Elements/elements.h

@@ -12 +12 @@
 
 double Paterson(double temperature);
+double Arrhenius(double temperature,double depth,double n);
 int    GetVerticesCoordinates(double* xyz,  Node** nodes, int numvertices);
 int    GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum);

issm/trunk/src/m/classes/model.m

@@ -139 +139 @@
                  rheology_B=NaN;
                  rheology_n=NaN;
+                 rheology_law=0;
 
                  %Geometrical parameters
@@ -732 +733 @@
                          md.hydro_kn=0;
 
+                         %Rheology law: what is the temperature dependence of B with T
+                         %available: None, Paterson and Arrhenius
+                         md.rheology_law=PatersonEnum;
 
                          %i/o:

issm/trunk/src/m/enum/EnumToString.m

@@ -376 +376 @@
 case OptionLogicalEnum(), string='OptionLogical'; return
 case OptionStructEnum(), string='OptionStruct'; return
+case RheologyLawEnum(), string='RheologyLaw'; return
+case PatersonEnum(), string='Paterson'; return
+case ArrheniusEnum(), string='Arrhenius'; return
 otherwise, error(['Enum ' num2str(enum)  ' not found']);
 

issm/trunk/src/m/enum/StringToEnum.m

@@ -374 +374 @@
                         elseif (strcmpi(name,'OptionLogical')), enum=OptionLogicalEnum(); return
                         elseif (strcmpi(name,'OptionStruct')), enum=OptionStructEnum(); return
+                        elseif (strcmpi(name,'RheologyLaw')), enum=RheologyLawEnum(); return
+                        elseif (strcmpi(name,'Paterson')), enum=PatersonEnum(); return
+                        elseif (strcmpi(name,'Arrhenius')), enum=ArrheniusEnum(); return
 else error(['Enum ' name  ' not found']);
 

issm/trunk/src/m/model/marshall.m

@@ -104 +104 @@
 WriteData(fid,md.rheology_B,'Mat','rheology_B');
 WriteData(fid,md.rheology_n,'Mat','rheology_n');
+WriteData(fid,md.rheology_law,'Integer','rheology_law');
 
 %Control methods