Changeset 9632

Timestamp:

09/06/11 16:36:00 (14 years ago)

Author:

seroussi

Message:

implemented thermal

Location:

issm/trunk

Files:

• src/c/EnumDefinitions/EnumDefinitions.h (modified) (4 diffs)
• src/c/modules/EnumToStringx/EnumToStringx.cpp (modified) (4 diffs)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp (modified) (1 diff)
• src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (modified) (1 diff)
• src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (modified) (2 diffs)
• src/c/modules/StringToEnumx/StringToEnumx.cpp (modified) (4 diffs)
• src/c/objects/Elements/Penta.cpp (modified) (18 diffs)
• src/c/objects/Loads/Pengrid.cpp (modified) (13 diffs)
• src/c/solvers/solver_thermal_nonlinear.cpp (modified) (4 diffs)
• src/m/classes/model/model.m (modified) (7 diffs)
• src/m/classes/thermal.m (added)
• src/m/enum/MinThermalConstraintsEnum.m (deleted)
• src/m/enum/SpctemperatureEnum.m (deleted)
• src/m/enum/StabilizeConstraintsEnum.m (deleted)
• src/m/enum/ThermalMaxiterEnum.m (added)
• src/m/enum/ThermalPenaltyFactorEnum.m (added)
• src/m/enum/ThermalPenaltyLockEnum.m (added)
• src/m/enum/ThermalPenaltyThresholdEnum.m (added)
• src/m/enum/ThermalRequestedOutputsEnum.m (added)
• src/m/enum/ThermalSpctemperatureEnum.m (added)
• src/m/enum/ThermalStabilizationEnum.m (added)
• src/m/model/collapse.m (modified) (1 diff)
• src/m/model/display/displaythermal.m (deleted)
• src/m/model/extrude.m (modified) (1 diff)
• src/m/model/ismodelselfconsistent.m (modified) (3 diffs)
• src/m/model/isresultconsistent.m (modified) (2 diffs)
• src/m/model/modelextract.m (modified) (1 diff)
• src/m/solvers/solver_thermal_nonlinear.m (modified) (2 diffs)
• src/m/utils/BC/SetIceSheetBC.m (modified) (1 diff)
• src/m/utils/BC/SetIceShelfBC.m (modified) (1 diff)
• src/m/utils/BC/SetMarineIceSheetBC.m (modified) (1 diff)
• test/Miscellaneous/GJM_test1/SquareShelf.par (modified) (1 diff)
• test/Miscellaneous/connectivity/Square.par (modified) (1 diff)
• test/NightlyRun/test1110.m (modified) (1 diff)
• test/NightlyRun/test1208.m (modified) (1 diff)
• test/NightlyRun/test1301.m (modified) (1 diff)
• test/NightlyRun/test1302.m (modified) (1 diff)
• test/NightlyRun/test1303.m (modified) (1 diff)
• test/NightlyRun/test1304.m (modified) (1 diff)
• test/NightlyRun/test263.m (modified) (1 diff)
• test/NightlyRun/test264.m (modified) (1 diff)
• test/NightlyRun/test265.m (modified) (1 diff)
• test/NightlyRun/test266.m (modified) (1 diff)
• test/NightlyRun/test529.m (modified) (1 diff)
• test/NightlyRun/test530.m (modified) (1 diff)
• test/NightlyRun/test531.m (modified) (1 diff)
• test/NightlyRun/test532.m (modified) (1 diff)
• test/Par/79North.par (modified) (1 diff)
• test/Par/ISMIPF.par (modified) (2 diffs)
• test/Par/RoundSheetShelf.par (modified) (1 diff)
• test/Par/SquareEISMINT.par (modified) (1 diff)
• test/Par/SquareSheetConstrained.par (modified) (1 diff)
• test/Par/SquareSheetShelf.par (modified) (1 diff)
• test/Par/SquareShelf.par (modified) (1 diff)
• test/Par/SquareShelfConstrained.par (modified) (1 diff)
• test/Par/SquareThermal.par (modified) (1 diff)

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

@@ -34 +34 @@
         SettingsLowmemEnum,
         SettingsIoGatherEnum,
+        ThermalSpctemperatureEnum,
+        ThermalPenaltyThresholdEnum,
+        ThermalPenaltyLockEnum,
+        ThermalMaxiterEnum,
+        ThermalStabilizationEnum,
+        ThermalPenaltyFactorEnum,
+        ThermalRequestedOutputsEnum,
         MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
         TimesteppingTimeStepEnum,
@@ -234 +241 @@
         SegmentOnIceShelfEnum,
         ShelfDampeningEnum,
-        StabilizeConstraintsEnum,
         SurfaceAreaEnum,
         SurfaceEnum,
@@ -350 +356 @@
         MaxNonlinearIterationsEnum,
         MinMechanicalConstraintsEnum,
-        MinThermalConstraintsEnum,
         NumberOfElementsEnum,
         NumberOfVerticesEnum,
@@ -425 +430 @@
         ZEnum,
         SpcthicknessEnum,
-        SpctemperatureEnum,
         PenaltyLockEnum,
         SpcvxEnum,

issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -38 +38 @@
                 case SettingsLowmemEnum : return "SettingsLowmem";
                 case SettingsIoGatherEnum : return "SettingsIoGather";
+                case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
+                case ThermalPenaltyThresholdEnum : return "ThermalPenaltyThreshold";
+                case ThermalPenaltyLockEnum : return "ThermalPenaltyLock";
+                case ThermalMaxiterEnum : return "ThermalMaxiter";
+                case ThermalStabilizationEnum : return "ThermalStabilization";
+                case ThermalPenaltyFactorEnum : return "ThermalPenaltyFactor";
+                case ThermalRequestedOutputsEnum : return "ThermalRequestedOutputs";
                 case MiscellaneousNameEnum : return "MiscellaneousName";
                 case TimesteppingTimeStepEnum : return "TimesteppingTimeStep";
@@ -201 +208 @@
                 case SegmentOnIceShelfEnum : return "SegmentOnIceShelf";
                 case ShelfDampeningEnum : return "ShelfDampening";
-                case StabilizeConstraintsEnum : return "StabilizeConstraints";
                 case SurfaceAreaEnum : return "SurfaceArea";
                 case SurfaceEnum : return "Surface";
@@ -301 +307 @@
                 case MaxNonlinearIterationsEnum : return "MaxNonlinearIterations";
                 case MinMechanicalConstraintsEnum : return "MinMechanicalConstraints";
-                case MinThermalConstraintsEnum : return "MinThermalConstraints";
                 case NumberOfElementsEnum : return "NumberOfElements";
                 case NumberOfVerticesEnum : return "NumberOfVertices";
@@ -368 +373 @@
                 case ZEnum : return "Z";
                 case SpcthicknessEnum : return "Spcthickness";
-                case SpctemperatureEnum : return "Spctemperature";
                 case PenaltyLockEnum : return "PenaltyLock";
                 case SpcvxEnum : return "Spcvx";

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

@@ -42 +42 @@
         parameters->AddObject(iomodel->CopyConstantObject(HydrostaticAdjustmentEnum));
         parameters->AddObject(iomodel->CopyConstantObject(PenaltyOffsetEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalPenaltyFactorEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SettingsLowmemEnum));
         parameters->AddObject(iomodel->CopyConstantObject(ConnectivityEnum));
@@ -51 +52 @@
         parameters->AddObject(iomodel->CopyConstantObject(ArtificialDiffusivityEnum));
         parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMeltingRateEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(MinThermalConstraintsEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalMaxiterEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalPenaltyThresholdEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalPenaltyLockEnum));
         parameters->AddObject(iomodel->CopyConstantObject(MinMechanicalConstraintsEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(StabilizeConstraintsEnum));
         parameters->AddObject(iomodel->CopyConstantObject(StokesreconditioningEnum));
         parameters->AddObject(iomodel->CopyConstantObject(ShelfDampeningEnum));

issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp

@@ -45 +45 @@
         /*Fetch data: */
         double *spctemperature=NULL;
-        iomodel->FetchData(&spctemperature,NULL,NULL,SpctemperatureEnum);
+        iomodel->FetchData(&spctemperature,NULL,NULL,ThermalSpctemperatureEnum);
 
         /*Initialize counter*/

issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp

@@ -33 +33 @@
         /*Only 3d mesh supported*/
         if (dim==3){
-                IoModelToConstraintsx(constraints,iomodel,SpctemperatureEnum,ThermalAnalysisEnum);
+                IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum);
         }
 

issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp

@@ -36 +36 @@
 
         //create penalties for nodes: no node can have a temperature over the melting point
-        iomodel->FetchData(2,SpctemperatureEnum,ElementsEnum);
+        iomodel->FetchData(2,ThermalSpctemperatureEnum,ElementsEnum);
         CreateSingleNodeToElementConnectivity(iomodel);
 
@@ -43 +43 @@
                 /*keep only this partition's nodes:*/
                 if((iomodel->my_vertices[i]==1)){
-                        if (isnan(iomodel->Data(SpctemperatureEnum)[i])){ //No penalty applied on spc nodes!
+                        if (isnan(iomodel->Data(ThermalSpctemperatureEnum)[i])){ //No penalty applied on spc nodes!
                                 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,ThermalAnalysisEnum));
                         }
                 }
         }
-        iomodel->DeleteData(2,SpctemperatureEnum,ElementsEnum);
+        iomodel->DeleteData(2,ThermalSpctemperatureEnum,ElementsEnum);
 
         /*Assign output pointer: */

issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -36 +36 @@
         else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
         else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
+        else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
+        else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
+        else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
+        else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
+        else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
+        else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
+        else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum;
         else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
         else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
@@ -199 +206 @@
         else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
         else if (strcmp(name,"ShelfDampening")==0) return ShelfDampeningEnum;
-        else if (strcmp(name,"StabilizeConstraints")==0) return StabilizeConstraintsEnum;
         else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
         else if (strcmp(name,"Surface")==0) return SurfaceEnum;
@@ -299 +305 @@
         else if (strcmp(name,"MaxNonlinearIterations")==0) return MaxNonlinearIterationsEnum;
         else if (strcmp(name,"MinMechanicalConstraints")==0) return MinMechanicalConstraintsEnum;
-        else if (strcmp(name,"MinThermalConstraints")==0) return MinThermalConstraintsEnum;
         else if (strcmp(name,"NumberOfElements")==0) return NumberOfElementsEnum;
         else if (strcmp(name,"NumberOfVertices")==0) return NumberOfVerticesEnum;
@@ -366 +371 @@
         else if (strcmp(name,"Z")==0) return ZEnum;
         else if (strcmp(name,"Spcthickness")==0) return SpcthicknessEnum;
-        else if (strcmp(name,"Spctemperature")==0) return SpctemperatureEnum;
         else if (strcmp(name,"PenaltyLock")==0) return PenaltyLockEnum;
         else if (strcmp(name,"Spcvx")==0) return SpcvxEnum;

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

@@ -1838 +1838 @@
 
         /*Intermediaries */
-        int        artdiff;
+        int        stabilization;
         int        i,j,ig,found=0;
         double     Jdet,u,v,w,um,vm,wm;
@@ -1852 +1852 @@
         double     B_conduct[3][numdof];
         double     B_advec[3][numdof];
-        double     B_artdiff[2][numdof];
+        double     B_stab[2][numdof];
         double     Bprime_advec[3][numdof];
         double     L[numdof];
         double     dbasis[3][6];
         double     D_scalar_conduct,D_scalar_advec;
-        double     D_scalar_trans,D_scalar_artdiff;
+        double     D_scalar_trans,D_scalar_stab;
         double     D[3][3];
         double     K[2][2]={0.0};
@@ -1875 +1875 @@
         thermalconductivity=matpar->GetThermalConductivity();
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->parameters->FindParam(&artdiff,ArtificialDiffusivityEnum);
+        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
         Input* pressure_input=inputs->GetInput(PressureEnum);      _assert_(pressure_input);
         Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);      _assert_(enthalpy_input);
@@ -1884 +1884 @@
         Input* vym_input=inputs->GetInput(VyMeshEnum);             _assert_(vym_input);
         Input* vzm_input=inputs->GetInput(VzMeshEnum);             _assert_(vzm_input);
-        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
+        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -1953 +1953 @@
                 /*Artifficial diffusivity*/
 
-                if(artdiff==1){
+                if(stabilization==1){
                         /*Build K: */
-                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
-                        if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
-                        K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
-                        K[1][0]=D_scalar_artdiff*fabs(u)*fabs(v);K[1][1]=D_scalar_artdiff*pow(v,2);
-
-                        GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
-
-                        TripleMultiply(&B_artdiff[0][0],2,numdof,1,
+                        D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
+                        if(dt) D_scalar_stab=D_scalar_stab*dt;
+                        K[0][0]=D_scalar_stab*pow(u,2);       K[0][1]=D_scalar_stab*fabs(u)*fabs(v);
+                        K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2);
+
+                        GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 
+
+                        TripleMultiply(&B_stab[0][0],2,numdof,1,
                                                 &K[0][0],2,2,0,
-                                                &B_artdiff[0][0],2,numdof,0,
+                                                &B_stab[0][0],2,numdof,0,
                                                 &Ke->values[0],1);
                 }
-                else if(artdiff==2){
+                else if(stabilization==2){
 
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
@@ -2115 +2115 @@
 
         /*Intermediaries */
-        int        artdiff;
+        int        stabilization;
         int        i,j,ig,found=0;
         double     Jdet,u,v,w,um,vm,wm;
@@ -2127 +2127 @@
         double     B_conduct[3][numdof];
         double     B_advec[3][numdof];
-        double     B_artdiff[2][numdof];
+        double     B_stab[2][numdof];
         double     Bprime_advec[3][numdof];
         double     L[numdof];
         double     dbasis[3][6];
         double     D_scalar_conduct,D_scalar_advec;
-        double     D_scalar_trans,D_scalar_artdiff;
+        double     D_scalar_trans,D_scalar_stab;
         double     D[3][3];
         double     K[2][2]={0.0};
@@ -2149 +2149 @@
         thermalconductivity=matpar->GetThermalConductivity();
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->parameters->FindParam(&artdiff,ArtificialDiffusivityEnum);
+        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
         Input* vx_input=inputs->GetInput(VxEnum);      _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum);      _assert_(vy_input);
@@ -2156 +2156 @@
         Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input);
         Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input);
-        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
+        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -2221 +2221 @@
                 /*Artifficial diffusivity*/
 
-                if(artdiff==1){
+                if(stabilization==1){
                         /*Build K: */
-                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
-                        if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
-                        K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
-                        K[1][0]=D_scalar_artdiff*fabs(u)*fabs(v);K[1][1]=D_scalar_artdiff*pow(v,2);
-
-                        GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
-
-                        TripleMultiply(&B_artdiff[0][0],2,numdof,1,
+                        D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
+                        if(dt) D_scalar_stab=D_scalar_stab*dt;
+                        K[0][0]=D_scalar_stab*pow(u,2);       K[0][1]=D_scalar_stab*fabs(u)*fabs(v);
+                        K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2);
+
+                        GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 
+
+                        TripleMultiply(&B_stab[0][0],2,numdof,1,
                                                 &K[0][0],2,2,0,
-                                                &B_artdiff[0][0],2,numdof,0,
+                                                &B_stab[0][0],2,numdof,0,
                                                 &Ke->values[0],1);
                 }
-                else if(artdiff==2){
+                else if(stabilization==2){
 
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
@@ -3287 +3287 @@
         /*Intermediaries*/
         int    i,j,ig,found=0;
-        int    friction_type,artdiff;
+        int    friction_type,stabilization;
         double Jdet,phi,dt;
         double rho_ice,heatcapacity;
@@ -3311 +3311 @@
         thermalconductivity=matpar->GetThermalConductivity();
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->parameters->FindParam(&artdiff,ArtificialDiffusivityEnum);
+        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -3317 +3317 @@
         Input* temperature_input=NULL;
         if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
-        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
+        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -3344 +3344 @@
                 }
 
-                if(artdiff==2){
+                if(stabilization==2){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
@@ -3554 +3554 @@
         /*Intermediaries*/
         int    i,j,ig,found=0;
-        int    friction_type,artdiff;
+        int    friction_type,stabilization;
         double Jdet,phi,dt;
         double rho_ice,heatcapacity;
@@ -3578 +3578 @@
         thermalconductivity=matpar->GetThermalConductivity();
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->parameters->FindParam(&artdiff,ArtificialDiffusivityEnum);
+        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -3584 +3584 @@
         Input* temperature_input=NULL;
         if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
-        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
+        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -3611 +3611 @@
                 }
 
-                if(artdiff==2){
+                if(stabilization==2){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 

issm/trunk/src/c/objects/Loads/Pengrid.cpp

@@ -476 +476 @@
         int    unstable=0;
         int    reset_penalties=0;
-        int    stabilize_constraints;
+        int    penalty_lock;
 
         /*recover pointers: */
@@ -493 +493 @@
 
         //Recover our data:
-        parameters->FindParam(&stabilize_constraints,StabilizeConstraintsEnum);
+        parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum);
         
         //Compute pressure melting point
@@ -514 +514 @@
         else{
                 unstable=1;
-                if(stabilize_constraints)zigzag_counter++;
+                if(penalty_lock)zigzag_counter++;
         }
 
         /*If penalty keeps zigzagging more than 5 times: */
-        if(stabilize_constraints){
-                if(zigzag_counter>stabilize_constraints){
+        if(penalty_lock){
+                if(zigzag_counter>penalty_lock){
                         unstable=0;
                         active=1;
@@ -566 +566 @@
         const int numdof=NUMVERTICES*NDOF1;
         double pressure,temperature,t_pmp;
-        double penalty_offset;
+        double penalty_factor;
 
         Penta* penta=(Penta*)element;
@@ -577 +577 @@
         penta->GetParameterValue(&pressure,node,PressureEnum);
         penta->GetParameterValue(&temperature,node,TemperatureEnum);
-        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
+        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
         
         /*Compute pressure melting point*/
@@ -584 +584 @@
         /*Add penalty load*/
         if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
-                Ke->values[0]=kmax*pow((double)10,penalty_offset);
+                Ke->values[0]=kmax*pow((double)10,penalty_factor);
         }
 
@@ -595 +595 @@
 
         const int numdof=NUMVERTICES*NDOF1;
-        double    penalty_offset;
+        double    penalty_factor;
 
         /*Initialize Element matrix and return if necessary*/
@@ -602 +602 @@
 
         /*recover parameters: */
-        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
-
-        Ke->values[0]=kmax*pow((double)10,penalty_offset);
+        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
+
+        Ke->values[0]=kmax*pow((double)10,penalty_factor);
 
         /*Clean up and return*/
@@ -618 +618 @@
         double melting_offset;
         double t_pmp;
-        double dt,penalty_offset;
+        double dt,penalty_factor;
 
         /*recover pointers: */
@@ -632 +632 @@
         inputs->GetParameterValue(&melting_offset,MeltingOffsetEnum);
         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
+        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
 
         /*Compute pressure melting point*/
@@ -640 +640 @@
           This time, the penalty must have the same value as the one used for the thermal computation
           so that the corresponding melting can be computed correctly
-          In the thermal computation, we used kmax=melting_offset, and the same penalty_offset*/
+          In the thermal computation, we used kmax=melting_offset, and the same penalty_factor*/
         if (temperature<t_pmp){ //%no melting
                 pe->values[0]=0;
         }
         else{
-                if (dt) pe->values[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
-                else    pe->values[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
+                if (dt) pe->values[0]=melting_offset*pow((double)10,penalty_factor)*(temperature-t_pmp)/dt;
+                else    pe->values[0]=melting_offset*pow((double)10,penalty_factor)*(temperature-t_pmp);
         }
 
@@ -659 +659 @@
         double pressure;
         double t_pmp;
-        double penalty_offset;
+        double penalty_factor;
 
         Penta* penta=(Penta*)element;
@@ -669 +669 @@
         /*Retrieve all inputs and parameters*/
         penta->GetParameterValue(&pressure,node,PressureEnum);
-        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
+        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
 
         /*Compute pressure melting point*/
         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
 
-        pe->values[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
+        pe->values[0]=kmax*pow((double)10,penalty_factor)*t_pmp;
 
         /*Clean up and return*/

issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp

@@ -27 +27 @@
         int num_unstable_constraints;
         int count;
-        int min_thermal_constraints;
+        int thermal_penalty_threshold;
+        int thermal_maxiter;
         bool reset_penalties;
 
@@ -39 +40 @@
         kflag=1; pflag=1;
         femmodel->parameters->FindParam(&lowmem,SettingsLowmemEnum);
-        femmodel->parameters->FindParam(&min_thermal_constraints,MinThermalConstraintsEnum);
+        femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
 
         count=1;
@@ -65 +67 @@
                 if (!converged){
                         _printf_(VerboseConvergence(),"%s%i\n","   #unstable constraints = ",num_unstable_constraints);
-                        if (num_unstable_constraints <= min_thermal_constraints)converged=true;
+                        if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
                 }
                 count++;
@@ -72 +74 @@
                 
                 if(converged)break;
+                if(count>=thermal_maxiter){
+                        _printf_(true,"   maximum number of iterations (%i) exceeded\n",thermal_maxiter); 
+                        break;
+                }
         }
 

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

@@ -23 +23 @@
                  settings  = modelfield('default',0,'marshall',true);
                  radaroverlay = modelfield('default',0,'marshall',false);
+                 thermal   = modelfield('default',0,'marshall',true);
                  miscellaneous = modelfield('default',0,'marshall',true);
                  timestepping = modelfield('default',0,'marshall',true);
@@ -111 +112 @@
                  mixed_layer_capacity       = modelfield('default',0,'marshall',true,'format','Double');
                  thermal_exchange_velocity  = modelfield('default',0,'marshall',true,'format','Double');
-                 min_thermal_constraints    = modelfield('default',0,'marshall',true,'format','Integer');
                  min_mechanical_constraints = modelfield('default',0,'marshall',true,'format','Integer');
-                 stabilize_constraints      = modelfield('default',0,'marshall',true,'format','Integer');
 
                  %Physical parameters
@@ -140 +139 @@
                  spcvy          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);
                  spcvz          = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);
-                 spctemperature = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);
                  spcthickness   = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);
                  diagnostic_ref = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);
@@ -409 +407 @@
                          if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end
                          if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end
+                         if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end
+                         if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end
+                         if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end
+                         if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end
+                         if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end
+                         if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end
                          if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end
                          if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end
@@ -530 +534 @@
                          md.settings=settings;
                          md.radaroverlay=radaroverlay;
+                         md.thermal=thermal;
                          md.miscellaneous=miscellaneous;
                          md.timestepping=timestepping;
@@ -621 +626 @@
                          %a few constraints remain unstable. For thermal computation, this
                          %parameter is often used.
-                         md.min_thermal_constraints=0;
                          md.min_mechanical_constraints=0;
 
@@ -709 +713 @@
                                  if(strcmp(index1.subs,'diagnostic')), displaydiagnostic(md);return; end
                                  if(strcmp(index1.subs,'prognostic')), displayprognostic(md);return; end
-                                 if(strcmp(index1.subs,'thermal')), displaythermal(md);return; end
                                  if(strcmp(index1.subs,'transient')), displaytransient(md);return; end
                                  if(strcmp(index1.subs,'control')), displaycontrol(md);return; end

issm/trunk/src/m/model/collapse.m

@@ -60 +60 @@
 md.spcvz=project2d(md,md.spcvz,md.numlayers);
 md.spcthickness=project2d(md,md.spcthickness,md.numlayers);
-md.spctemperature=project2d(md,md.spctemperature,md.numlayers);
+md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.numlayers);
 
 %Extrusion of Neumann BC

issm/trunk/src/m/model/extrude.m

@@ -175 +175 @@
 md.spcvy=project3d(md,'vector',md.spcvy,'type','node');
 md.spcvz=project3d(md,'vector',md.spcvz,'type','node');
-md.spctemperature=project3d(md,'vector',md.spctemperature,'type','node','layer',md.numlayers,'padding',NaN);
+md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.numlayers,'padding',NaN);
 md.spcthickness=project3d(md,'vector',md.spcthickness,'type','node');
 md.diagnostic_ref=project3d(md,'vector',md.diagnostic_ref,'type','node');

issm/trunk/src/m/model/ismodelselfconsistent.m

@@ -166 +166 @@
         message('model not consistent: artificial_diffusivity should be a scalar (0 or 1 or 2)');
 end
+if ~ismember(md.thermal.stabilization,[0 1 2]),
+        message('model not consistent: thermal.stabilization should be a scalar (0 or 1 or 2)');
+end
 if ~ismember(md.prognostic_DG,[0 1]),
         message('model not consistent: prognostic_DG should be a scalar (1 or 0)');
@@ -488 +491 @@
 
                         %CHECK THAT SPCTEMPERATURE IS AN APPROPRIATE FORCING
-                        fields={'spctemperature'};
+                        fields={'thermal.spctemperature'};
                         checkforcing(md,fields);
 
                         %CHECK THAT WE ARE NOT FULLY CONSTRAINED
-                        if ~any(~isnan(md.spctemperature))
+                        if ~any(~isnan(md.thermal.spctemperature))
                                 message(['model not consistent: model ' md.miscellaneous.name ' is totally constrained for temperature, no need to solve!']);
                         end
@@ -518 +521 @@
 
                                 %CHECK SPCTEMPERATURE that are not NaN are >0.
-                                if find(any(md.spctemperature(find(~isnan(md.spctemperature(1:md.numberofnodes,:))))<=0)),
+                                if find(any(md.thermal.spctemperature(find(~isnan(md.thermal.spctemperature(1:md.numberofnodes,:))))<=0)),
                                         message(['model not consistent: model ' md.miscellaneous.name ' is constrained with negative or nil temperatures!']);
                                 end

issm/trunk/src/m/model/isresultconsistent.m

@@ -147 +147 @@
                 %check melting (<=0 via penalties)
                 if any(abs(md.results.ThermalAnalysis(iter).melting(md.numberofnodes2d+1:end))>tolerance)
-                        disp(['''thermal'' result not consistent: increase penalty_offset (negative melting)']);
+                        disp(['''thermal'' result not consistent: increase thermal.penalty_factor (negative melting)']);
                         bool=0; return; 
                 end
@@ -197 +197 @@
                 if (md.dim==3),
                         if any(abs(md.results.TransientAnalysis(iter).melting(md.numberofnodes2d+1:end))>tolerance)
-                                disp(['''thermal'' result not consistent: increase penalty_offset (negative melting)']);
+                                disp(['''thermal'' result not consistent: increase therma.penalty_factor (negative melting)']);
                                 bool=0; return; 
                         end

issm/trunk/src/m/model/modelextract.m

@@ -208 +208 @@
                 md2.spcvz(nodestoflag2)=0;
         end
-        if ~isnan(md1.spctemperature),
-                md2.spctemperature(nodestoflag2,1)=1;
+        if ~isnan(md1.thermal.spctemperature),
+                md2.thermal.spctemperature(nodestoflag2,1)=1;
         end
 

issm/trunk/src/m/solvers/solver_thermal_nonlinear.m

@@ -11 +11 @@
         %Get parameters
         configuration_type=femmodel.parameters.ConfigurationType;
+        thermal_maxiter=femmodel.parameters.ThermalMaxiter;
 
         issmprintf(VerboseConvergence(),'\n%s',['   starting direct shooting method']);
@@ -36 +37 @@
                 if ~converged,
                         issmprintf(VerboseConvergence(),'%s%i','   #unstable constraints ',num_unstable_constraints);
-                        if num_unstable_constraints<=femmodel.parameters.MinThermalConstraints,
+                        if num_unstable_constraints<=femmodel.parameters.ThermalPenaltyThreshold,
                                 converged=1;
                         end
+                end
+                if(count>thermal_maxiter),
+                        issmprintf(true,'%s%i%s','      maximum number of iterations ',thermal_maxiter,' exceeded');
+                        break;
                 end
                 [femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);

issm/trunk/src/m/utils/BC/SetIceSheetBC.m

@@ -58 +58 @@
 
 if (length(md.temperature)==md.numberofnodes),
-        md.spctemperature=NaN*ones(md.numberofnodes,1);
-        pos=find(md.nodeonsurface); md.spctemperature(pos)=md.temperature(pos); %impose observed temperature on surface
+        md.thermal.spctemperature=NaN*ones(md.numberofnodes,1);
+        pos=find(md.nodeonsurface); md.thermal.spctemperature(pos)=md.temperature(pos); %impose observed temperature on surface
         if (length(md.basalforcings.geothermalflux)~=md.numberofnodes),
                 md.basalforcings.geothermalflux=50*10^-3*ones(md.numberofnodes,1); %50 mW/m^2

issm/trunk/src/m/utils/BC/SetIceShelfBC.m

@@ -90 +90 @@
 
 if (length(md.temperature)==md.numberofnodes),
-        md.spctemperature=NaN*ones(md.numberofnodes,1);
-        pos=find(md.nodeonsurface); md.spctemperature(pos)=md.temperature(pos); %impose observed temperature on surface
+        md.thermal.spctemperature=NaN*ones(md.numberofnodes,1);
+        pos=find(md.nodeonsurface); md.thermal.spctemperature(pos)=md.temperature(pos); %impose observed temperature on surface
         if (length(md.basalforcings.geothermalflux)~=md.numberofnodes),
                 md.basalforcings.geothermalflux=zeros(md.numberofnodes,1);

issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m

@@ -101 +101 @@
 
 if (length(md.temperature)==md.numberofnodes),
-        md.spctemperature=NaN*ones(md.numberofnodes,1);
-        pos=find(md.nodeonsurface); md.spctemperature(pos)=md.temperature(pos); %impose observed temperature on surface
+        md.thermal.spctemperature=NaN*ones(md.numberofnodes,1);
+        pos=find(md.nodeonsurface); md.thermal.spctemperature(pos)=md.temperature(pos); %impose observed temperature on surface
         if (length(md.basalforcings.geothermalflux)~=md.numberofnodes),
                 md.basalforcings.geothermalflux=zeros(md.numberofnodes,1);

issm/trunk/test/Miscellaneous/GJM_test1/SquareShelf.par

@@ -27 +27 @@
 md.viscosity_overshoot=0.3;
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 md.waitonlock=30;
 md.verbose=verbose(0);

issm/trunk/test/Miscellaneous/connectivity/Square.par

@@ -5 +5 @@
 md.ndt=md.dt*10; 
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 
 hmin=300;

issm/trunk/test/NightlyRun/test1110.m

@@ -36 +36 @@
         end
         pos=find(~md.nodeonbed);
-        md.spctemperature(pos)=255;
+        md.thermal.spctemperature(pos)=255;
 
         %Create MPCs to have periodic boundary conditions

issm/trunk/test/NightlyRun/test1208.m

@@ -23 +23 @@
 md.timestepping.final_time=50000;
 md.artificial_diffusivity=2;
+md.thermal.stabilization=2;
 
 %Now we can solve the problem 

issm/trunk/test/NightlyRun/test1301.m

@@ -16 +16 @@
 md.temperature=273.15*ones(md.numberofnodes,1);
 pos=find(md.nodeonsurface);
-md.spctemperature(pos)=md.temperature(pos);
+md.thermal.spctemperature(pos)=md.temperature(pos);
 md.rheology_B=paterson(md.temperature);
 

issm/trunk/test/NightlyRun/test1302.m

@@ -12 +12 @@
 
 %Thermal boundary conditions
-pos1=find(md.elementonbed);     md.spctemperature(md.elements(pos1,1:3))=10;
-pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6))=0;
+pos1=find(md.elementonbed);     md.thermal.spctemperature(md.elements(pos1,1:3))=10;
+pos2=find(md.elementonsurface); md.thermal.spctemperature(md.elements(pos2,4:6))=0;
 md.vz=0.1*ones(md.numberofnodes,1);
 md.vel=sqrt( md.vx.^2+ md.vy.^2+ md.vz.^2);

issm/trunk/test/NightlyRun/test1303.m

@@ -11 +11 @@
 md=extrude(md,11,2);
 md=setelementstype(md,'Pattyn','all');
-pos1=find(md.elementonbed);     md.spctemperature(md.elements(pos1,1:3))=10;
-pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6))=0;
+pos1=find(md.elementonbed);     md.thermal.spctemperature(md.elements(pos1,1:3))=10;
+pos2=find(md.elementonsurface); md.thermal.spctemperature(md.elements(pos2,4:6))=0;
 md.pressure=zeros(md.numberofnodes,1);
 

issm/trunk/test/NightlyRun/test1304.m

@@ -12 +12 @@
 md=setelementstype(md,'Pattyn','all');
 
-pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6))=0;
+pos2=find(md.elementonsurface); md.thermal.spctemperature(md.elements(pos2,4:6))=0;
 md.pressure=zeros(md.numberofnodes,1);
 md.basalforcings.geothermalflux(:)=0.1; %100mW/m^2

issm/trunk/test/NightlyRun/test263.m

@@ -5 +5 @@
 md=setelementstype(md,'macayeal','all');
 md.cluster=none;
-md.spctemperature=[md.spctemperature, md.spctemperature+5, md.spctemperature+10, md.spctemperature+15; 1.5 2.5 3.5 4];
+md.thermal.spctemperature=[md.thermal.spctemperature, md.thermal.spctemperature+5, md.thermal.spctemperature+10, md.thermal.spctemperature+15; 1.5 2.5 3.5 4];
 md.timestepping.time_step=1;
 md.timestepping.final_time=4;

issm/trunk/test/NightlyRun/test264.m

@@ -5 +5 @@
 md=setelementstype(md,'macayeal','all');
 md.cluster=generic('name',oshostname(),'np',3);
-md.spctemperature=[md.spctemperature, md.spctemperature+5, md.spctemperature+10, md.spctemperature+15; 1.5 2.5 3.5 4];
+md.thermal.spctemperature=[md.thermal.spctemperature, md.thermal.spctemperature+5, md.thermal.spctemperature+10, md.thermal.spctemperature+15; 1.5 2.5 3.5 4];
 md.timestepping.time_step=1;
 md.timestepping.final_time=4;

issm/trunk/test/NightlyRun/test265.m

@@ -5 +5 @@
 md=setelementstype(md,'pattyn','all');
 md.cluster=none;
-md.spctemperature=[md.spctemperature, md.spctemperature+5; 1 2];
+md.thermal.spctemperature=[md.thermal.spctemperature, md.thermal.spctemperature+5; 1 2];
 md.timestepping.time_step=0.5;
 md.timestepping.final_time=2;

issm/trunk/test/NightlyRun/test266.m

@@ -5 +5 @@
 md=setelementstype(md,'pattyn','all');
 md.cluster=generic('name',oshostname(),'np',3);
-md.spctemperature=[md.spctemperature, md.spctemperature+5; 1 2];
+md.thermal.spctemperature=[md.thermal.spctemperature, md.thermal.spctemperature+5; 1 2];
 md.timestepping.time_step=0.5;
 md.timestepping.final_time=2;

issm/trunk/test/NightlyRun/test529.m

@@ -4 +4 @@
 md=extrude(md,3,1);
 md=setelementstype(md,'pattyn','all');
-md.artificial_diffusivity=2;
+md.thermal.stabilization=2;
 md=solve(md,ThermalSolutionEnum);
 

issm/trunk/test/NightlyRun/test530.m

@@ -4 +4 @@
 md=extrude(md,3,1);
 md=setelementstype(md,'pattyn','all');
-md.artificial_diffusivity=2;
+md.thermal.stabilization=2;
 md.cluster=generic('name',oshostname(),'np',3);
 md=solve(md,ThermalSolutionEnum);

issm/trunk/test/NightlyRun/test531.m

@@ -4 +4 @@
 md=extrude(md,3,1);
 md=setelementstype(md,'pattyn','all');
-md.artificial_diffusivity=2;
+md.thermal.stabilization=2;
 md.timestepping.time_step=0;
-md.min_thermal_constraints=40;
+md.thermal.penalty_threshold=40;
 md=solve(md,ThermalSolutionEnum);
 

issm/trunk/test/NightlyRun/test532.m

@@ -4 +4 @@
 md=extrude(md,3,1);
 md=setelementstype(md,'pattyn','all');
-md.artificial_diffusivity=2;
+md.thermal.stabilization=2;
 md.cluster=generic('name',oshostname(),'np',3);
 md.timestepping.time_step=0;
-md.min_thermal_constraints=40;
+md.thermal.penalty_threshold=40;
 md=solve(md,ThermalSolutionEnum);
 

issm/trunk/test/Par/79North.par

@@ -32 +32 @@
 md.viscosity_overshoot=0.3;
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 md.verbose=verbose(0);
 md.waitonlock=30;

issm/trunk/test/Par/ISMIPF.par

@@ -31 +31 @@
 md.spcthickness=NaN*ones(md.numberofnodes,1);
 md.spcthickness(pos)=md.thickness(pos);
-md.spctemperature=255*ones(md.numberofnodes,1);
+md.thermal.spctemperature=255*ones(md.numberofnodes,1);
 md.basalforcings.geothermalflux=0.4*ones(md.numberofnodes,1);
 
@@ -41 +41 @@
 md.timestepping.final_time=10;
 md.artificial_diffusivity=1;
-md.min_thermal_constraints=10^5;
+md.thermal.stabilization=1;
+md.thermal.penalty_threshold=10^5;

issm/trunk/test/Par/RoundSheetShelf.par

@@ -59 +59 @@
 md.viscosity_overshoot=0.0;
 md.artificial_diffusivity=1;
+md.thermal.stabilisation=1;
 md.verbose=verbose(0);
 md.waitonlock=30;

issm/trunk/test/Par/SquareEISMINT.par

@@ -34 +34 @@
 md.spcthickness=NaN*ones(md.numberofnodes,1);
 md.spcthickness(pos)=500;
+md.artificial_diffusivity=0; %Better result with no artificial diffusivity
+md.thermal.statbilization=0; 
 md.timestepping.final_time=500;
 md.timestepping.time_step=1;
-md.artificial_diffusivity=0; %Better result with no artificial diffusivity

issm/trunk/test/Par/SquareSheetConstrained.par

@@ -33 +33 @@
 md.viscosity_overshoot=0.0;
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 md.verbose=verbose(0);
 md.waitonlock=30;

issm/trunk/test/Par/SquareSheetShelf.par

@@ -40 +40 @@
 md.viscosity_overshoot=0.0;
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 md.verbose=verbose(0);
 md.waitonlock=30;

issm/trunk/test/Par/SquareShelf.par

@@ -33 +33 @@
 md.viscosity_overshoot=0.3;
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 md.waitonlock=30;
 md.verbose=verbose(0);

issm/trunk/test/Par/SquareShelfConstrained.par

@@ -37 +37 @@
 md.viscosity_overshoot=0.0;
 md.artificial_diffusivity=1;
+md.thermal.stabilization=1;
 md.verbose=verbose(0);
 md.waitonlock=30;

issm/trunk/test/Par/SquareThermal.par

@@ -39 +39 @@
 
 disp('      boundary conditions for thermal model');
-md.spctemperature(:)=md.temperature;
+md.thermal.spctemperature(:)=md.temperature;
 md.basalforcings.geothermalflux=zeros(md.numberofnodes,1); 
 pos=find(md.elementonicesheet);md.basalforcings.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2