Changeset 15079

Timestamp:

05/22/13 10:22:19 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: base convergence on vectors rather than inputs (more accurate and makes more sense)

File:

• issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/analyses/steadystate_core.cpp

@@ -16 +16 @@
 
 /*Local prototypes*/
-bool steadystateconvergence(FemModel* femmodel);
+bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol);
 
 void steadystate_core(FemModel* femmodel){
@@ -22 +22 @@
         /*intermediary: */
         int step; 
+        Vector<IssmDouble>* ug     = NULL;
+        Vector<IssmDouble>* ug_old = NULL;
+        Vector<IssmDouble>* tg     = NULL;
+        Vector<IssmDouble>* tg_old = NULL;
 
         /*parameters: */
-        bool save_results,isenthalpy;
-        int  maxiter;
-        int  numoutputs         = 0;
-        int  *requested_outputs = NULL;
+        bool        save_results,isenthalpy;
+        int         maxiter;
+        IssmDouble  reltol;
+        int         numoutputs        = 0;
+        int        *requested_outputs = NULL;
 
         /* recover parameters:*/
@@ -34 +39 @@
         femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
         femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+        femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum);
         femmodel->parameters->SetParam(false,SaveResultsEnum);
         if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);
@@ -46 +52 @@
                 if(isenthalpy==0){
                         thermal_core(femmodel);
+                        femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
+                        GetSolutionFromInputsx(&tg,femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
                 }
                 else{
                         enthalpy_core(femmodel);
+                        GetSolutionFromInputsx(&tg,femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
                 }
                 #else
@@ -56 +65 @@
                 if(VerboseSolution()) _pprintLine_("   computing new velocity");
                 diagnostic_core(femmodel);
+                GetSolutionFromInputsx(&ug,femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 
-                if (step>1){
-                        if(VerboseSolution()) _pprintLine_("   checking velocity, temperature and pressure convergence");
-                        if(steadystateconvergence(femmodel)) break;
+                if(step>1){
+                        if(VerboseSolution()) _pprintLine_("   checking steadystate convergence");
+                        if(steadystateconvergence(tg,tg_old,ug,ug_old,reltol)) break;
                 }
                 if(step>maxiter){
@@ -66 +76 @@
                 }
 
-                if(VerboseSolution()) _pprintLine_("   saving velocity, temperature and pressure to check for convergence at next step");
-                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxPicardEnum);
-                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyPicardEnum);
-                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzPicardEnum);
-                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressurePicardEnum);
-                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum);
-
-                //increase counter
+                /*update results and increase counter*/
+                delete tg_old;tg_old=tg;
+                delete ug_old;ug_old=ug;
                 step++;
         }
@@ -92 +97 @@
 
         /*Free ressources:*/
+        delete tg;
+        delete ug;
         xDelete<int>(requested_outputs);
 }
-bool steadystateconvergence(FemModel* femmodel){
+bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol){
 
-        /*output: */
-        bool converged=false;
-        bool velocity_converged=false;
-        bool temperature_converged=false;
+        /*Output*/
+        bool converged = true;
 
-        /*intermediary: */
-        int velocityenums[8]={VxEnum,VxPicardEnum,VyEnum,VyPicardEnum,VzEnum,VzPicardEnum,PressureEnum,PressurePicardEnum}; //pairs of enums (new and old) on which to carry out the converence tests
-        int temperatureenums[2]={TemperatureEnum,TemperatureOldEnum};
-        int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum
-        IssmDouble convergencecriterionvalue[1]; //value of criterion to be respected
+        /*Intermediary*/
+        Vector<IssmDouble>* dug    = NULL;
+        Vector<IssmDouble>* dtg    = NULL;
+        IssmDouble          ndt,nt;
+        IssmDouble          ndu,nu;
 
-        /*retrieve parameters: */
-        femmodel->parameters->FindParam(&convergencecriterionvalue[0],SteadystateReltolEnum);
+        /*compute norm(du)/norm(u)*/
+        dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
+        ndu=dug->Norm(NORM_TWO); nu=ug_old->Norm(NORM_TWO);
+        if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+        if((ndu/nu)<reltol){
+                if(VerboseConvergence()) _pprintLine_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" < "<<reltol*100<<" %");
+        }
+        else{ 
+                if(VerboseConvergence()) _pprintLine_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" > "<<reltol*100<<" %");
+                converged=false;
+        }
 
-        /*figure out convergence at the input level, because we don't have the solution vectors!: */
-        velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&velocityenums[0],8,&convergencecriterion[0],&convergencecriterionvalue[0],1);
-        temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&temperatureenums[0],2,&convergencecriterion[0],&convergencecriterionvalue[0],1);
+        /*compute norm(dt)/norm(t)*/
+        dtg=tg_old->Duplicate(); tg_old->Copy(dtg); dtg->AYPX(tg,-1.0);
+        ndt=dtg->Norm(NORM_TWO); nt=tg_old->Norm(NORM_TWO);
+        if (xIsNan<IssmDouble>(ndt) || xIsNan<IssmDouble>(nt)) _error_("convergence criterion is NaN!");
+        if((ndt/nt)<reltol){
+                if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" < "<<reltol*100<<" %");
+        }
+        else{ 
+                if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" > "<<reltol*100<<" %");
+                converged=false;
+        }
 
-        if(velocity_converged && temperature_converged) converged=true;
-
-        /*return: */
+        /*clean up and return*/
+        delete dtg;
+        delete dug;
         return converged;
 }