Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 9631) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 9632) @@ -34,4 +34,11 @@ 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,5 +241,4 @@ SegmentOnIceShelfEnum, ShelfDampeningEnum, - StabilizeConstraintsEnum, SurfaceAreaEnum, SurfaceEnum, @@ -350,5 +356,4 @@ MaxNonlinearIterationsEnum, MinMechanicalConstraintsEnum, - MinThermalConstraintsEnum, NumberOfElementsEnum, NumberOfVerticesEnum, @@ -425,5 +430,4 @@ ZEnum, SpcthicknessEnum, - SpctemperatureEnum, PenaltyLockEnum, SpcvxEnum, Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 9631) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 9632) @@ -38,4 +38,11 @@ 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,5 +208,4 @@ case SegmentOnIceShelfEnum : return "SegmentOnIceShelf"; case ShelfDampeningEnum : return "ShelfDampening"; - case StabilizeConstraintsEnum : return "StabilizeConstraints"; case SurfaceAreaEnum : return "SurfaceArea"; case SurfaceEnum : return "Surface"; @@ -301,5 +307,4 @@ case MaxNonlinearIterationsEnum : return "MaxNonlinearIterations"; case MinMechanicalConstraintsEnum : return "MinMechanicalConstraints"; - case MinThermalConstraintsEnum : return "MinThermalConstraints"; case NumberOfElementsEnum : return "NumberOfElements"; case NumberOfVerticesEnum : return "NumberOfVertices"; @@ -368,5 +373,4 @@ case ZEnum : return "Z"; case SpcthicknessEnum : return "Spcthickness"; - case SpctemperatureEnum : return "Spctemperature"; case PenaltyLockEnum : return "PenaltyLock"; case SpcvxEnum : return "Spcvx"; Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 9631) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 9632) @@ -42,4 +42,5 @@ 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,7 +52,9 @@ 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)); Index: /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp (revision 9631) +++ /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp (revision 9632) @@ -45,5 +45,5 @@ /*Fetch data: */ double *spctemperature=NULL; - iomodel->FetchData(&spctemperature,NULL,NULL,SpctemperatureEnum); + iomodel->FetchData(&spctemperature,NULL,NULL,ThermalSpctemperatureEnum); /*Initialize counter*/ Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 9631) +++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 9632) @@ -33,5 +33,5 @@ /*Only 3d mesh supported*/ if (dim==3){ - IoModelToConstraintsx(constraints,iomodel,SpctemperatureEnum,ThermalAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum); } Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 9631) +++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 9632) @@ -36,5 +36,5 @@ //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,10 +43,10 @@ /*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: */ Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 9631) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 9632) @@ -36,4 +36,11 @@ 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,5 +206,4 @@ 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,5 +305,4 @@ 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,5 +371,4 @@ 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; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 9631) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 9632) @@ -1838,5 +1838,5 @@ /*Intermediaries */ - int artdiff; + int stabilization; int i,j,ig,found=0; double Jdet,u,v,w,um,vm,wm; @@ -1852,10 +1852,10 @@ 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,5 +1875,5 @@ 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,5 +1884,5 @@ 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,19 +1953,19 @@ /*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,5 +2115,5 @@ /*Intermediaries */ - int artdiff; + int stabilization; int i,j,ig,found=0; double Jdet,u,v,w,um,vm,wm; @@ -2127,10 +2127,10 @@ 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,5 +2149,5 @@ 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,5 +2156,5 @@ 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,19 +2221,19 @@ /*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,5 +3287,5 @@ /*Intermediaries*/ int i,j,ig,found=0; - int friction_type,artdiff; + int friction_type,stabilization; double Jdet,phi,dt; double rho_ice,heatcapacity; @@ -3311,5 +3311,5 @@ 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,5 +3317,5 @@ 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,5 +3344,5 @@ } - if(artdiff==2){ + if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); @@ -3554,5 +3554,5 @@ /*Intermediaries*/ int i,j,ig,found=0; - int friction_type,artdiff; + int friction_type,stabilization; double Jdet,phi,dt; double rho_ice,heatcapacity; @@ -3578,5 +3578,5 @@ 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,5 +3584,5 @@ 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,5 +3611,5 @@ } - if(artdiff==2){ + if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); Index: /issm/trunk/src/c/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 9631) +++ /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 9632) @@ -476,5 +476,5 @@ int unstable=0; int reset_penalties=0; - int stabilize_constraints; + int penalty_lock; /*recover pointers: */ @@ -493,5 +493,5 @@ //Recover our data: - parameters->FindParam(&stabilize_constraints,StabilizeConstraintsEnum); + parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum); //Compute pressure melting point @@ -514,10 +514,10 @@ 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,5 +566,5 @@ const int numdof=NUMVERTICES*NDOF1; double pressure,temperature,t_pmp; - double penalty_offset; + double penalty_factor; Penta* penta=(Penta*)element; @@ -577,5 +577,5 @@ 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,5 +584,5 @@ /*Add penalty load*/ if (temperaturevalues[0]=kmax*pow((double)10,penalty_offset); + Ke->values[0]=kmax*pow((double)10,penalty_factor); } @@ -595,5 +595,5 @@ const int numdof=NUMVERTICES*NDOF1; - double penalty_offset; + double penalty_factor; /*Initialize Element matrix and return if necessary*/ @@ -602,7 +602,7 @@ /*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,5 +618,5 @@ double melting_offset; double t_pmp; - double dt,penalty_offset; + double dt,penalty_factor; /*recover pointers: */ @@ -632,5 +632,5 @@ inputs->GetParameterValue(&melting_offset,MeltingOffsetEnum); parameters->FindParam(&dt,TimesteppingTimeStepEnum); - parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); + parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); /*Compute pressure melting point*/ @@ -640,11 +640,11 @@ 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 (temperaturevalues[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,5 +659,5 @@ double pressure; double t_pmp; - double penalty_offset; + double penalty_factor; Penta* penta=(Penta*)element; @@ -669,10 +669,10 @@ /*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*/ Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 9631) +++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 9632) @@ -27,5 +27,6 @@ int num_unstable_constraints; int count; - int min_thermal_constraints; + int thermal_penalty_threshold; + int thermal_maxiter; bool reset_penalties; @@ -39,6 +40,7 @@ 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,5 +67,5 @@ 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,4 +74,8 @@ if(converged)break; + if(count>=thermal_maxiter){ + _printf_(true," maximum number of iterations (%i) exceeded\n",thermal_maxiter); + break; + } } Index: /issm/trunk/src/m/classes/model/model.m =================================================================== --- /issm/trunk/src/m/classes/model/model.m (revision 9631) +++ /issm/trunk/src/m/classes/model/model.m (revision 9632) @@ -23,4 +23,5 @@ 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,7 +112,5 @@ 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,5 +139,4 @@ 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,4 +407,10 @@ 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,4 +534,5 @@ md.settings=settings; md.radaroverlay=radaroverlay; + md.thermal=thermal; md.miscellaneous=miscellaneous; md.timestepping=timestepping; @@ -621,5 +626,4 @@ %a few constraints remain unstable. For thermal computation, this %parameter is often used. - md.min_thermal_constraints=0; md.min_mechanical_constraints=0; @@ -709,5 +713,4 @@ 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 Index: /issm/trunk/src/m/classes/thermal.m =================================================================== --- /issm/trunk/src/m/classes/thermal.m (revision 9632) +++ /issm/trunk/src/m/classes/thermal.m (revision 9632) @@ -0,0 +1,65 @@ +%THERMAL class definition +% +% Usage: +% thermal=thermal(); + +classdef thermal + properties (SetAccess=public) + spctemperature = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1); + penalty_threshold = modelfield('default',0,'marshall',true,'format','Integer'); + stabilization = modelfield('default',0,'marshall',true,'format','Integer'); + maxiter = modelfield('default',0,'marshall',true,'format','Integer'); + penalty_lock = modelfield('default',0,'marshall',true,'format','Integer'); + penalty_factor = modelfield('default',0,'marshall',true,'format','Double'); + requested_outputs = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3); + end + methods + function obj = thermal(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + in=varargin{1}; + if (isa(in,'numeric') & in==0), + % requesting templates do nothing + else + error('constructor not supported'); + end + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + %first, use the defaults provided by the properties definition above. + fieldnames=fields(obj); + for i=1:length(fieldnames), + fieldname=fieldnames{i}; + obj.(fieldname)=obj.(fieldname).default; + end + + %Number of unstable constraints acceptable + obj.penalty_threshold=0; + + %Type of stabilization used + obj.stabilization=1; + + %Maximum number of iterations + obj.maxiter=100; + + %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor + obj.penalty_factor=3; + end % }}} + function disp(obj) % {{{ + disp(sprintf(' Thermal solution parameters:')); + + disp(sprintf('\n parameters:')); + fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); + fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)'); + + disp(sprintf('\n boundary conditions:')); + fielddisplay(obj,'spctemperature','constraints flag list (first column) and values (second column)'); + + end % }}} + end +end Index: sm/trunk/src/m/enum/MinThermalConstraintsEnum.m =================================================================== --- /issm/trunk/src/m/enum/MinThermalConstraintsEnum.m (revision 9631) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=MinThermalConstraintsEnum() -%MINTHERMALCONSTRAINTSENUM - Enum of MinThermalConstraints -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh -% Please read src/c/EnumDefinitions/README for more information -% -% Usage: -% macro=MinThermalConstraintsEnum() - -macro=StringToEnum('MinThermalConstraints'); Index: sm/trunk/src/m/enum/SpctemperatureEnum.m =================================================================== --- /issm/trunk/src/m/enum/SpctemperatureEnum.m (revision 9631) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SpctemperatureEnum() -%SPCTEMPERATUREENUM - Enum of Spctemperature -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh -% Please read src/c/EnumDefinitions/README for more information -% -% Usage: -% macro=SpctemperatureEnum() - -macro=StringToEnum('Spctemperature'); Index: sm/trunk/src/m/enum/StabilizeConstraintsEnum.m =================================================================== --- /issm/trunk/src/m/enum/StabilizeConstraintsEnum.m (revision 9631) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=StabilizeConstraintsEnum() -%STABILIZECONSTRAINTSENUM - Enum of StabilizeConstraints -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh -% Please read src/c/EnumDefinitions/README for more information -% -% Usage: -% macro=StabilizeConstraintsEnum() - -macro=StringToEnum('StabilizeConstraints'); Index: /issm/trunk/src/m/enum/ThermalMaxiterEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalMaxiterEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalMaxiterEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalMaxiterEnum() +%THERMALMAXITERENUM - Enum of ThermalMaxiter +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalMaxiterEnum() + +macro=StringToEnum('ThermalMaxiter'); Index: /issm/trunk/src/m/enum/ThermalPenaltyFactorEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalPenaltyFactorEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalPenaltyFactorEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalPenaltyFactorEnum() +%THERMALPENALTYFACTORENUM - Enum of ThermalPenaltyFactor +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalPenaltyFactorEnum() + +macro=StringToEnum('ThermalPenaltyFactor'); Index: /issm/trunk/src/m/enum/ThermalPenaltyLockEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalPenaltyLockEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalPenaltyLockEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalPenaltyLockEnum() +%THERMALPENALTYLOCKENUM - Enum of ThermalPenaltyLock +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalPenaltyLockEnum() + +macro=StringToEnum('ThermalPenaltyLock'); Index: /issm/trunk/src/m/enum/ThermalPenaltyThresholdEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalPenaltyThresholdEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalPenaltyThresholdEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalPenaltyThresholdEnum() +%THERMALPENALTYTHRESHOLDENUM - Enum of ThermalPenaltyThreshold +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalPenaltyThresholdEnum() + +macro=StringToEnum('ThermalPenaltyThreshold'); Index: /issm/trunk/src/m/enum/ThermalRequestedOutputsEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalRequestedOutputsEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalRequestedOutputsEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalRequestedOutputsEnum() +%THERMALREQUESTEDOUTPUTSENUM - Enum of ThermalRequestedOutputs +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalRequestedOutputsEnum() + +macro=StringToEnum('ThermalRequestedOutputs'); Index: /issm/trunk/src/m/enum/ThermalSpctemperatureEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalSpctemperatureEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalSpctemperatureEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalSpctemperatureEnum() +%THERMALSPCTEMPERATUREENUM - Enum of ThermalSpctemperature +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalSpctemperatureEnum() + +macro=StringToEnum('ThermalSpctemperature'); Index: /issm/trunk/src/m/enum/ThermalStabilizationEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThermalStabilizationEnum.m (revision 9632) +++ /issm/trunk/src/m/enum/ThermalStabilizationEnum.m (revision 9632) @@ -0,0 +1,11 @@ +function macro=ThermalStabilizationEnum() +%THERMALSTABILIZATIONENUM - Enum of ThermalStabilization +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalStabilizationEnum() + +macro=StringToEnum('ThermalStabilization'); Index: /issm/trunk/src/m/model/collapse.m =================================================================== --- /issm/trunk/src/m/model/collapse.m (revision 9631) +++ /issm/trunk/src/m/model/collapse.m (revision 9632) @@ -60,5 +60,5 @@ 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 Index: sm/trunk/src/m/model/display/displaythermal.m =================================================================== --- /issm/trunk/src/m/model/display/displaythermal.m (revision 9631) +++ (revision ) @@ -1,22 +1,0 @@ -function displaythermal(md) -%DISPLAYTHERMAL - display solution parameters applied to thermal model -% -% To avoid clobbering display.m with every field from model md, -% we create this routine that displays the solution parameters -% from model md, only if requested. -% -% Usage: -% displaythermal(md) - -disp(sprintf(' Thermal solution parameters:')); - -disp(sprintf('\n parameters:')); -fielddisplay(md,'stabilize_constraints','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); -fielddisplay(md,'min_thermal_constraints','threshold to declare convergence of thermal solution (default is 0)'); - -disp(sprintf('\n boundary conditions:')); -fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)'); - -disp(sprintf('\n transient:')); -fielddisplay(md,'dt','time step [yr]'); -fielddisplay(md,'ndt','time span [yr]'); Index: /issm/trunk/src/m/model/extrude.m =================================================================== --- /issm/trunk/src/m/model/extrude.m (revision 9631) +++ /issm/trunk/src/m/model/extrude.m (revision 9632) @@ -175,5 +175,5 @@ 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'); Index: /issm/trunk/src/m/model/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 9631) +++ /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 9632) @@ -166,4 +166,7 @@ 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,9 +491,9 @@ %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,5 +521,5 @@ %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 Index: /issm/trunk/src/m/model/isresultconsistent.m =================================================================== --- /issm/trunk/src/m/model/isresultconsistent.m (revision 9631) +++ /issm/trunk/src/m/model/isresultconsistent.m (revision 9632) @@ -147,5 +147,5 @@ %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,5 +197,5 @@ 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 Index: /issm/trunk/src/m/model/modelextract.m =================================================================== --- /issm/trunk/src/m/model/modelextract.m (revision 9631) +++ /issm/trunk/src/m/model/modelextract.m (revision 9632) @@ -208,6 +208,6 @@ 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 Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m =================================================================== --- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m (revision 9631) +++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m (revision 9632) @@ -11,4 +11,5 @@ %Get parameters configuration_type=femmodel.parameters.ConfigurationType; + thermal_maxiter=femmodel.parameters.ThermalMaxiter; issmprintf(VerboseConvergence(),'\n%s',[' starting direct shooting method']); @@ -36,7 +37,11 @@ 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); Index: /issm/trunk/src/m/utils/BC/SetIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetIceSheetBC.m (revision 9631) +++ /issm/trunk/src/m/utils/BC/SetIceSheetBC.m (revision 9632) @@ -58,6 +58,6 @@ 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 Index: /issm/trunk/src/m/utils/BC/SetIceShelfBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetIceShelfBC.m (revision 9631) +++ /issm/trunk/src/m/utils/BC/SetIceShelfBC.m (revision 9632) @@ -90,6 +90,6 @@ 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); Index: /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 9631) +++ /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 9632) @@ -101,6 +101,6 @@ 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); Index: /issm/trunk/test/Miscellaneous/GJM_test1/SquareShelf.par =================================================================== --- /issm/trunk/test/Miscellaneous/GJM_test1/SquareShelf.par (revision 9631) +++ /issm/trunk/test/Miscellaneous/GJM_test1/SquareShelf.par (revision 9632) @@ -27,4 +27,5 @@ md.viscosity_overshoot=0.3; md.artificial_diffusivity=1; +md.thermal.stabilization=1; md.waitonlock=30; md.verbose=verbose(0); Index: /issm/trunk/test/Miscellaneous/connectivity/Square.par =================================================================== --- /issm/trunk/test/Miscellaneous/connectivity/Square.par (revision 9631) +++ /issm/trunk/test/Miscellaneous/connectivity/Square.par (revision 9632) @@ -5,4 +5,5 @@ md.ndt=md.dt*10; md.artificial_diffusivity=1; +md.thermal.stabilization=1; hmin=300; Index: /issm/trunk/test/NightlyRun/test1110.m =================================================================== --- /issm/trunk/test/NightlyRun/test1110.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test1110.m (revision 9632) @@ -36,5 +36,5 @@ end pos=find(~md.nodeonbed); - md.spctemperature(pos)=255; + md.thermal.spctemperature(pos)=255; %Create MPCs to have periodic boundary conditions Index: /issm/trunk/test/NightlyRun/test1208.m =================================================================== --- /issm/trunk/test/NightlyRun/test1208.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test1208.m (revision 9632) @@ -23,4 +23,5 @@ md.timestepping.final_time=50000; md.artificial_diffusivity=2; +md.thermal.stabilization=2; %Now we can solve the problem Index: /issm/trunk/test/NightlyRun/test1301.m =================================================================== --- /issm/trunk/test/NightlyRun/test1301.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test1301.m (revision 9632) @@ -16,5 +16,5 @@ 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); Index: /issm/trunk/test/NightlyRun/test1302.m =================================================================== --- /issm/trunk/test/NightlyRun/test1302.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test1302.m (revision 9632) @@ -12,6 +12,6 @@ %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); Index: /issm/trunk/test/NightlyRun/test1303.m =================================================================== --- /issm/trunk/test/NightlyRun/test1303.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test1303.m (revision 9632) @@ -11,6 +11,6 @@ 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); Index: /issm/trunk/test/NightlyRun/test1304.m =================================================================== --- /issm/trunk/test/NightlyRun/test1304.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test1304.m (revision 9632) @@ -12,5 +12,5 @@ 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 Index: /issm/trunk/test/NightlyRun/test263.m =================================================================== --- /issm/trunk/test/NightlyRun/test263.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test263.m (revision 9632) @@ -5,5 +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; Index: /issm/trunk/test/NightlyRun/test264.m =================================================================== --- /issm/trunk/test/NightlyRun/test264.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test264.m (revision 9632) @@ -5,5 +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; Index: /issm/trunk/test/NightlyRun/test265.m =================================================================== --- /issm/trunk/test/NightlyRun/test265.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test265.m (revision 9632) @@ -5,5 +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; Index: /issm/trunk/test/NightlyRun/test266.m =================================================================== --- /issm/trunk/test/NightlyRun/test266.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test266.m (revision 9632) @@ -5,5 +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; Index: /issm/trunk/test/NightlyRun/test529.m =================================================================== --- /issm/trunk/test/NightlyRun/test529.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test529.m (revision 9632) @@ -4,5 +4,5 @@ md=extrude(md,3,1); md=setelementstype(md,'pattyn','all'); -md.artificial_diffusivity=2; +md.thermal.stabilization=2; md=solve(md,ThermalSolutionEnum); Index: /issm/trunk/test/NightlyRun/test530.m =================================================================== --- /issm/trunk/test/NightlyRun/test530.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test530.m (revision 9632) @@ -4,5 +4,5 @@ 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); Index: /issm/trunk/test/NightlyRun/test531.m =================================================================== --- /issm/trunk/test/NightlyRun/test531.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test531.m (revision 9632) @@ -4,7 +4,7 @@ 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); Index: /issm/trunk/test/NightlyRun/test532.m =================================================================== --- /issm/trunk/test/NightlyRun/test532.m (revision 9631) +++ /issm/trunk/test/NightlyRun/test532.m (revision 9632) @@ -4,8 +4,8 @@ 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); Index: /issm/trunk/test/Par/79North.par =================================================================== --- /issm/trunk/test/Par/79North.par (revision 9631) +++ /issm/trunk/test/Par/79North.par (revision 9632) @@ -32,4 +32,5 @@ md.viscosity_overshoot=0.3; md.artificial_diffusivity=1; +md.thermal.stabilization=1; md.verbose=verbose(0); md.waitonlock=30; Index: /issm/trunk/test/Par/ISMIPF.par =================================================================== --- /issm/trunk/test/Par/ISMIPF.par (revision 9631) +++ /issm/trunk/test/Par/ISMIPF.par (revision 9632) @@ -31,5 +31,5 @@ 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,3 +41,4 @@ 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; Index: /issm/trunk/test/Par/RoundSheetShelf.par =================================================================== --- /issm/trunk/test/Par/RoundSheetShelf.par (revision 9631) +++ /issm/trunk/test/Par/RoundSheetShelf.par (revision 9632) @@ -59,4 +59,5 @@ md.viscosity_overshoot=0.0; md.artificial_diffusivity=1; +md.thermal.stabilisation=1; md.verbose=verbose(0); md.waitonlock=30; Index: /issm/trunk/test/Par/SquareEISMINT.par =================================================================== --- /issm/trunk/test/Par/SquareEISMINT.par (revision 9631) +++ /issm/trunk/test/Par/SquareEISMINT.par (revision 9632) @@ -34,5 +34,6 @@ 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 Index: /issm/trunk/test/Par/SquareSheetConstrained.par =================================================================== --- /issm/trunk/test/Par/SquareSheetConstrained.par (revision 9631) +++ /issm/trunk/test/Par/SquareSheetConstrained.par (revision 9632) @@ -33,4 +33,5 @@ md.viscosity_overshoot=0.0; md.artificial_diffusivity=1; +md.thermal.stabilization=1; md.verbose=verbose(0); md.waitonlock=30; Index: /issm/trunk/test/Par/SquareSheetShelf.par =================================================================== --- /issm/trunk/test/Par/SquareSheetShelf.par (revision 9631) +++ /issm/trunk/test/Par/SquareSheetShelf.par (revision 9632) @@ -40,4 +40,5 @@ md.viscosity_overshoot=0.0; md.artificial_diffusivity=1; +md.thermal.stabilization=1; md.verbose=verbose(0); md.waitonlock=30; Index: /issm/trunk/test/Par/SquareShelf.par =================================================================== --- /issm/trunk/test/Par/SquareShelf.par (revision 9631) +++ /issm/trunk/test/Par/SquareShelf.par (revision 9632) @@ -33,4 +33,5 @@ md.viscosity_overshoot=0.3; md.artificial_diffusivity=1; +md.thermal.stabilization=1; md.waitonlock=30; md.verbose=verbose(0); Index: /issm/trunk/test/Par/SquareShelfConstrained.par =================================================================== --- /issm/trunk/test/Par/SquareShelfConstrained.par (revision 9631) +++ /issm/trunk/test/Par/SquareShelfConstrained.par (revision 9632) @@ -37,4 +37,5 @@ md.viscosity_overshoot=0.0; md.artificial_diffusivity=1; +md.thermal.stabilization=1; md.verbose=verbose(0); md.waitonlock=30; Index: /issm/trunk/test/Par/SquareThermal.par =================================================================== --- /issm/trunk/test/Par/SquareThermal.par (revision 9631) +++ /issm/trunk/test/Par/SquareThermal.par (revision 9632) @@ -39,5 +39,5 @@ 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