Changeset 11299

Timestamp:

02/02/12 07:45:03 (13 years ago)

Author:

seroussi

Message:

added non linearity in enthalpy solution

Location:

issm/trunk-jpl/src

Files:

• c/objects/Elements/Penta.cpp (modified) (5 diffs)
• c/objects/Materials/Matpar.cpp (modified) (1 diff)
• c/solutions/enthalpy_core.cpp (modified) (1 diff)
• c/solutions/enthalpy_core_step.cpp (modified) (1 diff)
• m/solutions/enthalpy_core.m (modified) (2 diffs)
• m/solutions/enthalpy_core_step.m (added)
• m/solvers/solver_nonlinear.m (modified) (3 diffs)

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

@@ -2050 +2050 @@
                                 name==InversionVzObsEnum ||
                                 name==TemperatureEnum ||
+                                name==EnthalpyEnum ||
+                                name==EnthalpyPicardEnum ||
                                 name==WaterfractionEnum||
                                 name==FrictionCoefficientEnum ||
@@ -3241 +3243 @@
         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
         Input* pressure_input=inputs->GetInput(PressureEnum);      _assert_(pressure_input);
-        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);      _assert_(enthalpy_input);
+        Input* enthalpy_input=inputs->GetInput(EnthalpyPicardEnum);_assert_(enthalpy_input); //for this iteration of the step
         Input* vx_input=inputs->GetInput(VxEnum);                  _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum);                  _assert_(vy_input);
@@ -4239 +4241 @@
         const int    numdof=NDOF1*NUMVERTICES;
 
-        bool   converged;
+        bool   converged=false;
         int    i,rheology_law;
         double xyz_list[NUMVERTICES][3];
@@ -4266 +4268 @@
         Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
         
-
-//      this->inputs->GetInputValue(&converged,ConvergedEnum);
-//      if(converged){
+        this->inputs->GetInputValue(&converged,ConvergedEnum);
+        if(converged){
                 /*Convert enthalpy into temperature and water fraction*/
                 for(i=0;i<numdof;i++) matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
@@ -4300 +4301 @@
 
                 }
-//      }
-//      else{
-//              this->inputs->AddInput(new PentaP1Input(EnthalpyPicardEnum,values));
-//      }
+        }
+        else{
+                this->inputs->AddInput(new PentaP1Input(EnthalpyPicardEnum,values));
+        }
 
         /*Free ressources:*/

issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp ¶

@@ -420 +420 @@
         }
         else{
+                printf("-------------- file: Matpar.cpp line: %i\n",__LINE__); 
                 return 0.1*thermalconductivity/(rho_ice*heatcapacity);
         }

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

@@ -46 +46 @@
 
                 /*call enthalpy_core_step: */
-                femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
-                solver_linear(femmodel);
+                enthalpy_core_step(femmodel,i,time);
 
                 if(solution_type==EnthalpySolutionEnum && !control_analysis){

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

@@ -21 +21 @@
         solver_nonlinear(femmodel,modify_loads);
 
-        /*transfer water column thickness to old water column thickness: */
-        //InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,EnthalpyOldEnum);
+        /*transfer enthalpy to enthalpy picard for the next step: */
+        InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,EnthalpyPicardEnum);
 
 }

issm/trunk-jpl/src/m/solutions/enthalpy_core.m ¶

@@ -4 +4 @@
 %   Usage:
 %      femmodel=enthalpy_core(femmodel)
-
 
         %recover parameters common to all solutions
@@ -25 +24 @@
                 time=(i+1)*dt;
 
-                issmprintf(VerboseSolution,'\n%s',['   computing enthalpy']);
-                femmodel=SetCurrentConfiguration(femmodel,EnthalpyAnalysisEnum);
-                femmodel=solver_linear(femmodel);
+                femmodel=enthalpy_core_step(femmodel,i,time); 
 
                 if (solution_type==EnthalpySolutionEnum & ~control_analysis),

issm/trunk-jpl/src/m/solvers/solver_nonlinear.m ¶

@@ -22 +22 @@
 
         %Update the solution to make sure that vx and vxold are similar
+        [femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
         [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
 
@@ -38 +39 @@
                 ug= Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); 
 
-                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
-                [loads,constraints_converged,num_unstable_constraints] =ConstraintsState( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters);
-
-                issmprintf(VerboseConvergence(),'%s%i','      number of unstable constraints: ',num_unstable_constraints);
-
                 %Figure out if convergence have been reached
                 converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
@@ -48 +44 @@
                 %add convergence status into  status
                 [femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
+
+                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
+                [loads,constraints_converged,num_unstable_constraints] =ConstraintsState( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters);
+
+                issmprintf(VerboseConvergence(),'%s%i','      number of unstable constraints: ',num_unstable_constraints);
 
                 %rift convergence