Changeset 18553

Timestamp:

10/01/14 05:58:08 (10 years ago)

Author:

jbondzio

Message:

ADD: thermal.reltol for enthalpy convergence check

Location:

issm/trunk-jpl/src

Files:

• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• m/classes/thermal.m (modified) (6 diffs)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/ThermalReltolEnum.m (added)

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -274 +274 @@
         ThermalIsenthalpyEnum,
         ThermalIsdynamicbasalspcEnum,
+        ThermalReltolEnum,
         ThermalMaxiterEnum,
         ThermalPenaltyFactorEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -282 +282 @@
                 case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
                 case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
+                case ThermalReltolEnum : return "ThermalReltol";
                 case ThermalMaxiterEnum : return "ThermalMaxiter";
                 case ThermalPenaltyFactorEnum : return "ThermalPenaltyFactor";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -288 +288 @@
               else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
               else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
+              else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
               else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
               else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
@@ -382 +383 @@
               else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
               else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
-              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+              if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
+              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
               else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
               else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
@@ -505 +506 @@
               else if (strcmp(name,"Fill")==0) return FillEnum;
               else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
-              else if (strcmp(name,"Friction")==0) return FrictionEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"Internal")==0) return InternalEnum;
+              if (strcmp(name,"Friction")==0) return FrictionEnum;
+              else if (strcmp(name,"Internal")==0) return InternalEnum;
               else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
               else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
@@ -628 +629 @@
               else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
               else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
-              else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+              if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
+              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
               else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
               else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
@@ -751 +752 @@
               else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum;
               else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
-              else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
+              if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
+              else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
               else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
               else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;

issm/trunk-jpl/src/m/classes/thermal.m

@@ -9 +9 @@
                 penalty_threshold = 0;
                 stabilization     = 0;
+                reltol                          = 0;
                 maxiter           = 0;
                 penalty_lock      = 0;
@@ -31 +32 @@
 
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="',class(obj.stabilization),'" default="',convert2str(obj.stabilization),'">','     <section name="thermal" />','     <help> maximum number of non linear iterations </help>','  </parameter>');
+                                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="reltol" type="',class(obj.reltol),'" default="',convert2str(obj.reltol),'">','     <section name="steadystate" />','     <help> relative tolerance criterion [K] </help>','  </parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter" type="',class(obj.maxiter),'" default="',convert2str(obj.maxiter),'">','     <section name="thermal" />','     <help> stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','  </parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_lock" type="',class(obj.penalty_lock),'" default="',convert2str(obj.penalty_lock),'">','     <section name="thermal" />','     <help> threshold to declare convergence of thermal solution (default is 0)  </help>','  </parameter>');
@@ -65 +67 @@
                         obj.stabilization=1;
 
+                        %Relative tolerance for the enthalpy convergence
+                        obj.reltol=0.01;
+
                         %Maximum number of iterations
                         obj.maxiter=100;
@@ -94 +99 @@
                                 md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[0 1]);
                                 if(md.thermal.isenthalpy)
+                                        md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message','reltol must be larger than zero');
+                                        if isnan(md.stressbalance.reltol),
+                                                md = checkmessage(md,['for a steadystate computation, stressbalance.reltol (relative convergence criterion) must be defined!']);
+                                        end
                                         md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[1], 'message',['for enthalpy run thermal.isdynamicbasalspc should be 1']);
                                 end
@@ -105 +114 @@
                         fielddisplay(obj,'spctemperature','temperature constraints (NaN means no constraint) [K]');
                         fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
+                        fielddisplay(obj,'reltol','relative tolerance convergence criterion for enthalpy');
                         fielddisplay(obj,'maxiter','maximum number of non linear iterations');
                         fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
@@ -118 +128 @@
                         WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
                         WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
+                        WriteData(fid,'object',obj,'fieldname','reltol','format','Double');
                         WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
                         WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -274 +274 @@
 def ThermalIsenthalpyEnum(): return StringToEnum("ThermalIsenthalpy")[0]
 def ThermalIsdynamicbasalspcEnum(): return StringToEnum("ThermalIsdynamicbasalspc")[0]
+def ThermalReltolEnum(): return StringToEnum("ThermalReltol")[0]
 def ThermalMaxiterEnum(): return StringToEnum("ThermalMaxiter")[0]
 def ThermalPenaltyFactorEnum(): return StringToEnum("ThermalPenaltyFactor")[0]