Changeset 16745

Timestamp:

11/13/13 16:15:51 (11 years ago)

Author:

seroussi

Message:

CHG: working on Thermal and Enthalpy input update from solution

Location:

issm/trunk-jpl/src/c

Files:

• analyses/EnthalpyAnalysis.cpp (modified) (5 diffs)
• analyses/ThermalAnalysis.cpp (modified) (6 diffs)
• classes/Elements/Penta.cpp (modified) (2 diffs)

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

@@ -202 +202 @@
         IssmDouble* values        = xNew<IssmDouble>(numnodes);
         IssmDouble* pressure      = xNew<IssmDouble>(numnodes);
+        IssmDouble* surface       = xNew<IssmDouble>(numnodes);
+        IssmDouble* B             = xNew<IssmDouble>(numnodes);
         IssmDouble* temperature   = xNew<IssmDouble>(numnodes);
         IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
@@ -215 +217 @@
         /*Get all inputs and parameters*/
         element->GetInputValue(&converged,ConvergedEnum);
-        element->GetInputListOnNodes(pressure,PressureEnum);
+        element->GetInputListOnNodes(&pressure[0],PressureEnum);
         if(converged){
                 for(i=0;i<numnodes;i++){
@@ -229 +231 @@
                  * otherwise the rheology could be negative*/
                 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+                element->GetInputListOnNodes(&surface[0],SurfaceEnum);
                 switch(rheology_law){
                         case NoneEnum:
@@ -234 +237 @@
                                 break;
                         case PatersonEnum:
-                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-                                B_average=Paterson(T_average);
-                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+                                for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
                                 break;
-                        case ArrheniusEnum:{
-                                Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
-                                surface_input->GetInputAverage(&s_average);
+                        case ArrheniusEnum:
                                 element->GetVerticesCoordinates(&xyz_list);
-                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-                                //B_average=Arrhenius(T_average,
-                                                        //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),
-                                                        //element->GetMaticeParameter(MaterialsRheologyNEnum));
-                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
                                 break;
-                                }
                         case LliboutryDuvalEnum:
-                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-                                for(i=0;i<numnodes;i++) P_average+=pressure[i]/reCast<IssmDouble>(numnodes);
-                                B_average=LliboutryDuval(T_average,P_average,
-                                                        element->GetMaterialParameter(MaterialsRheologyNEnum),
-                                                        element->GetMaterialParameter(MaterialsBetaEnum),
-                                                        element->GetMaterialParameter(ConstantsReferencetemperatureEnum),
-                                                        element->GetMaterialParameter(MaterialsHeatcapacityEnum),
-                                                        element->GetMaterialParameter(MaterialsLatentheatEnum));
-                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
-                                break;
-                        default:
-                                _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
+                                //for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat()); for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]); 
+                                for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum); 
+                                break; 
+                        default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
                 }
         }
@@ -271 +260 @@
         xDelete<IssmDouble>(values);
         xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(surface);
+        xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(temperature);
         xDelete<IssmDouble>(waterfraction);

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

@@ -119 +119 @@
         bool        converged;
         int         i,rheology_law;
-        IssmDouble  B_average,s_average,T_average=0.;
         int        *doflist   = NULL;
         IssmDouble *xyz_list  = NULL;
@@ -130 +129 @@
         element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
         IssmDouble* values    = xNew<IssmDouble>(numnodes);
+        IssmDouble* surface   = xNew<IssmDouble>(numnodes);
+        IssmDouble* B         = xNew<IssmDouble>(numnodes);
 
         /*Use the dof list to index into the solution vector: */
@@ -144 +145 @@
         if(hack){
                 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
-                element->GetInputListOnVertices(pressure,PressureEnum);
+                element->GetInputListOnNodes(&pressure[0],PressureEnum);
                 for(i=0;i<numnodes;i++){
                         if(values[i]>element->TMeltingPoint(pressure[i]))     values[i]=element->TMeltingPoint(pressure[i]);
@@ -160 +161 @@
                  * otherwise the rheology could be negative*/
                 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+                element->GetInputListOnNodes(&surface[0],SurfaceEnum);
                 switch(rheology_law){
                         case NoneEnum:
@@ -165 +167 @@
                                 break;
                         case PatersonEnum:
-                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-                                B_average=Paterson(T_average);
-                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+                                for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
                                 break;
                         case ArrheniusEnum:{
-                                Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
-                                surface_input->GetInputAverage(&s_average);
                                 element->GetVerticesCoordinates(&xyz_list);
-                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-                                //B_average=Arrhenius(T_average,
-                                                        //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),
-                                                        //element->GetMaticeParameter(MaterialsRheologyNEnum));
-                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
                                 break;
                                 }
                         default:
                                 _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
-
                 }
         }
@@ -191 +186 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(surface);
+        xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(xyz_list);
         xDelete<int>(doflist);

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

@@ -4929 +4929 @@
         IssmDouble waterfraction[numdof];
         IssmDouble B[numdof];
-        IssmDouble B_average,s_average;
         int*       doflist=NULL;
 
@@ -4964 +4963 @@
                  * otherwise the rheology could be negative*/
                 this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
-                GetInputListOnVertices(&surface[0],PressureEnum);
+                GetInputListOnVertices(&surface[0],SurfaceEnum);
                 switch(rheology_law){
                         case NoneEnum: