Changeset 9013

Timestamp:

07/16/11 01:38:53 (14 years ago)

Author:

Eric.Larour

Message:

Added MaxSteadystateIterations, to control steady state a little bit better.

Location:

issm/trunk/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (1 diff)
• c/objects/IoModel.cpp (modified) (2 diffs)
• c/objects/IoModel.h (modified) (1 diff)
• c/solutions/steadystate_core.cpp (modified) (3 diffs)
• m/classes/model.m (modified) (2 diffs)
• m/enum/MaxSteadystateIterationsEnum.m (added)
• m/model/extrude.m (modified) (1 diff)
• m/model/partition/AreaAverageOntoPartition.m (modified) (2 diffs)
• m/model/partition/partitioner.m (modified) (2 diffs)
• m/model/tres.m (modified) (1 diff)

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

@@ -541 +541 @@
         NpartEnum,
         QmuMassFluxNumProfilesEnum,
-        PartEnum
+        PartEnum,
+        MaxSteadystateIterationsEnum
 };
 

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

@@ -483 +483 @@
                 case QmuMassFluxNumProfilesEnum : return "QmuMassFluxNumProfiles";
                 case PartEnum : return "Part";
+                case MaxSteadystateIterationsEnum : return "MaxSteadystateIterations";
                 default : return "unknown";
 

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

@@ -39 +39 @@
         parameters->AddObject(new DoubleParam(EpsAbsEnum,iomodel->eps_abs));
         parameters->AddObject(new IntParam(MaxNonlinearIterationsEnum,(IssmInt)iomodel->max_nonlinear_iterations));
+        parameters->AddObject(new IntParam(MaxSteadystateIterationsEnum,(IssmInt)iomodel->max_steadystate_iterations));
         parameters->AddObject(new DoubleParam(EpsvelEnum,iomodel->epsvel));
         parameters->AddObject(new DoubleParam(YtsEnum,iomodel->yts));

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

@@ -481 +481 @@
         else if (strcmp(name,"QmuMassFluxNumProfiles")==0) return QmuMassFluxNumProfilesEnum;
         else if (strcmp(name,"Part")==0) return PartEnum;
+        else if (strcmp(name,"MaxSteadystateIterations")==0) return MaxSteadystateIterationsEnum;
         else _error_("Enum %s not found",name);
 

issm/trunk/src/c/objects/IoModel.cpp

@@ -185 +185 @@
         IoModelFetchData(&this->eps_abs,iomodel_handle,EpsAbsEnum);
         IoModelFetchData(&this->max_nonlinear_iterations,iomodel_handle,MaxNonlinearIterationsEnum);
+        IoModelFetchData(&this->max_steadystate_iterations,iomodel_handle,MaxSteadystateIterationsEnum);
         IoModelFetchData(&this->dt,iomodel_handle,DtEnum);
         IoModelFetchData(&this->ndt,iomodel_handle,NdtEnum);
@@ -355 +356 @@
         this->eps_abs=0;
         this->max_nonlinear_iterations=0;
+        this->max_steadystate_iterations=0;
         this->dt=0;
         this->ndt=0;

issm/trunk/src/c/objects/IoModel.h

@@ -150 +150 @@
                 double  eps_abs;
                 double  max_nonlinear_iterations;
+                double  max_steadystate_iterations;
                 double  dt,ndt;
                 int     time_adapt;

issm/trunk/src/c/solutions/steadystate_core.cpp

@@ -20 +20 @@
         int dim;
         int solution_type;
-        int max_its;
+        int max_steadystate_iterations;
         bool control_analysis;
         
@@ -27 +27 @@
         femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-        femmodel->parameters->FindParam(&max_its,MaxNonlinearIterationsEnum);
+        femmodel->parameters->FindParam(&max_steadystate_iterations,MaxSteadystateIterationsEnum);
 
         /*intialize counters: */
@@ -44 +44 @@
                         if(steadystateconvergence(femmodel)) break;
                 }
-                if(step>max_its){
-                        _printf_(VerboseSolution(),"%s%i%s\n","   maximum number of iterations ",max_its," reached");
+                if(step>max_steadystate_iterations){
+                        _printf_(VerboseSolution(),"%s%i%s\n","   maximum number steadystate iterations ",max_steadystate_iterations," reached");
                         break;
                 }

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

@@ -193 +193 @@
                  eps_abs                  = {0,true,'Double'};
                  max_nonlinear_iterations = {0,true,'Double'};
+                 max_steadystate_iterations = {0,true,'Double'};
                  sparsity                 = {0,true,'Double'};
                  connectivity             = {0,true,'Integer'};
@@ -660 +661 @@
                          %maximum of non-linear iterations.
                          md.max_nonlinear_iterations=100;
+                         
+                         %maximum of steady state iterations
+                         md.max_steadystate_iterations=100;
 
                          %sparsity

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

@@ -224 +224 @@
 md.surface=project3d(md,'vector',md.surface,'type','node');
 md.thickness=project3d(md,'vector',md.thickness,'type','node');
+md.thickness_coeff=project3d(md,'vector',md.thickness_coeff,'type','node');
 md.bed=project3d(md,'vector',md.bed,'type','node');
 md.nodeonboundary=project3d(md,'vector',md.nodeonboundary,'type','node');

issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m

@@ -1 @@
-function partvector=AreaAverageOntoPartition(md,vector)
+function partvector=AreaAverageOntoPartition(md,vector,layer)
 %AREAAVERAGEONTOPARTITION  compute partition values for a certain vector expressed on the vertices of the mesh. Use area weighted average.
 %
 %   Usage: average=AreaAverageOntoPartition(md,vector)
+%           average=AreaAverageOntoPartition(md,vector,layer) %if in 3D, chose which layer is partitioned
 %
 
-%ok, first check that part is Matlab matlab indexed
+%some checks
+if md.dim==3,
+        if nargin~=3,
+                error('layer should be provided onto which Area Averaging occurs');
+        end
+        %save 3D model
+        md3d=md;
+        
+        md.elements=md.elements2d;
+        md.x=md.x2d;
+        md.y=md.y2d;
+        md.numberofnodes=md.numberofnodes2d;
+        md.numberofelements=md.numberofelements2d;
+        md.vwgt=[];
+        md.nodeconnectivity=[];
+
+        %run connectivity routine
+        md=adjacency(md);
+
+        %finally, project vector: 
+        vector=project2d(md3d,vector,layer);
+        md.part=project2d(md3d,md3d.part,layer);
+end
+
+%ok, first check that part is Matlab indexed
 part=md.part+1;
+
+%some check: 
+if md.npart~=max(part),
+        error('AreaAverageOntoPartition error message: ''npart'' should be equal to max(md.part)');
+end
 
 %initialize output
@@ -17 +47 @@
         partvector(i)=sum(weightedvector(pos))/sum(md.vwgt(pos));
 end
+
+%in 3D, restore 3D model:
+if md.dim==3,
+        md=md3d;
+end

issm/trunk/src/m/model/partition/partitioner.m

@@ -30 +30 @@
 npart=getfieldvalue(options,'npart');
 recomputeadjacency=getfieldvalue(options,'recomputeadjacency');
+
+if(md.dim==3),
+        %partitioning essentially happens in 2D. So partition in 2D, then 
+        %extrude the partition vector vertically. 
+        md3d=md; %save  for later
+        md.elements=md.elements2d;
+        md.x=md.x2d;
+        md.y=md.y2d;
+        md.numberofnodes=md.numberofnodes2d;
+        md.numberofelements=md.numberofelements2d;
+        md.vwgt=[];
+        md.nodeconnectivity=[];
+end
 
 %adjacency matrix if needed:
@@ -87 +100 @@
 end
 
+%extrude if we are in 3D:
+if md.dim==3,
+        md=md3d;
+        part=project3d(md,'vector',part','type','node');
+end
+
 md.part=part;

issm/trunk/src/m/model/tres.m

@@ -85 +85 @@
         md.temperature=PatchToVec(md.results.SteadystateSolution.Temperature);
         md.basal_melting_rate=PatchToVec(md.results.SteadystateSolution.BasalMeltingRate);
-        
+
         if md.control_analysis==1,
-                if control_type==md.control_type
-                        md.(EnumToModelField(control_type))=PatchToVec(md.results.DiagnosticSolution.(EnumToString(control_type)));
+                for control_type=md.control_type
+                        md.(EnumToModelField(control_type))=PatchToVec(md.results.SteadystateSolution.(EnumToString(control_type)));
                 end
         end