Changeset 21382

Timestamp:

11/16/16 19:51:51 (8 years ago)

Author:

felicity

Message:

NEW: ESTAR thermal model capability

Location:

issm/trunk-jpl/src/c/analyses

Files:

• EnthalpyAnalysis.cpp (modified) (9 diffs)
• ThermalAnalysis.cpp (modified) (7 diffs)

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

@@ -108 +108 @@
 
         bool dakota_analysis,ismovingfront,isenthalpy;
-        int frictionlaw,basalforcing_model;
+        int frictionlaw,basalforcing_model,materialstype;
         int FrictionCoupling;
         
@@ -134 +134 @@
         iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
         iomodel->FindConstant(&frictionlaw,"md.friction.law");
+        iomodel->FindConstant(&materialstype,"md.materials.type");
 
         iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
@@ -145 +146 @@
                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
         }
-        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
-        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
@@ -171 +170 @@
                         iomodel->FetchDataToInput(elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
                         break;
+        }
+
+        /*Rheology type*/
+        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
+        switch(materialstype){
+                case MatdamageiceEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        break;
+                case MatestarEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
+                        break;
+                case MaticeEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        break;
+                default:
+                        _error_("not supported");
         }
         
@@ -589 +605 @@
                 }
 
-                /*Artifficial diffusivity*/
+                /*Artificial diffusivity*/
                 if(stabilization==1){
                         element->ElementSizes(&hx,&hy,&hz);
@@ -1452 +1468 @@
         int         i,rheology_law;
         IssmDouble  B_average,s_average,T_average=0.,P_average=0.;
+        IssmDouble  n=3.0;
         int        *doflist   = NULL;
         IssmDouble *xyz_list  = NULL;
@@ -1477 +1494 @@
 
         /*Get all inputs and parameters*/
+        if(element->material->ObjectEnum()!=MatestarEnum) n=element->GetMaterialParameter(MaterialsRheologyNEnum);
         element->GetInputValue(&converged,ConvergedEnum);
         element->GetInputListOnNodes(&pressure[0],PressureEnum);
@@ -1506 +1524 @@
                                 break;
                         case CuffeyTemperateEnum:
-                                for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i], element->GetMaterialParameter(MaterialsRheologyNEnum));
+                                for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i],n);
                                 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
                                 break;
@@ -1515 +1533 @@
                         case ArrheniusEnum:
                                 element->GetVerticesCoordinates(&xyz_list);
-                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],n);
                                 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
                                 break;
                         case LliboutryDuvalEnum:
-                                for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],element->GetMaterialParameter(MaterialsRheologyNEnum),element->GetMaterialParameter(MaterialsBetaEnum),element->GetMaterialParameter(ConstantsReferencetemperatureEnum),element->GetMaterialParameter(MaterialsHeatcapacityEnum),element->GetMaterialParameter(MaterialsLatentheatEnum)); 
+                                for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],n,element->GetMaterialParameter(MaterialsBetaEnum),element->GetMaterialParameter(ConstantsReferencetemperatureEnum),element->GetMaterialParameter(MaterialsHeatcapacityEnum),element->GetMaterialParameter(MaterialsLatentheatEnum)); 
                                 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); 
                                 break; 

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

@@ -73 +73 @@
 void ThermalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
-        int frictionlaw,basalforcing_model;
+        int frictionlaw,basalforcing_model,materialstype;
         int FrictionCoupling;
         /*Now, is the model 3d? otherwise, do nothing: */
@@ -93 +93 @@
         iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
         iomodel->FindConstant(&frictionlaw,"md.friction.law");
+        iomodel->FindConstant(&materialstype,"md.materials.type");
 
         iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
@@ -106 +107 @@
         iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
         iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
-        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
-        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
@@ -116 +115 @@
         InputUpdateFromConstantx(elements,0.,VyMeshEnum);
         InputUpdateFromConstantx(elements,0.,VzMeshEnum);
+
+        /*Rheology type*/
+        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
+        switch(materialstype){
+                case MatdamageiceEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        break;
+                case MatestarEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
+                        break;
+                case MaticeEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        break;
+                default:
+                        _error_("not supported");
+        }
         if(ismovingfront){
                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
@@ -722 +738 @@
         int        *doflist   = NULL;
         IssmDouble *xyz_list  = NULL;
+        IssmDouble  n=3.0;
         bool        hack      = false;
 
@@ -756 +773 @@
 
         /*Get all inputs and parameters*/
+        if(element->material->ObjectEnum()!=MatestarEnum) n=element->GetMaterialParameter(MaterialsRheologyNEnum);
         element->GetInputValue(&converged,ConvergedEnum);
         if(converged){
@@ -782 +800 @@
                         case ArrheniusEnum:{
                                 element->GetVerticesCoordinates(&xyz_list);
-                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],n);
                                 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
                                 break;