Changeset 18207

Timestamp:

07/01/14 10:37:09 (11 years ago)

Author:

seroussi

Message:

BUG: fixing Uzawa pressure

Location:

issm/trunk-jpl/src/c

Files:

• analyses/UzawaPressureAnalysis.cpp (modified) (4 diffs)
• solutionsequences/solutionsequence_la.cpp (modified) (7 diffs)

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

@@ -92 +92 @@
         int          dim;
         IssmDouble   Jdet,r,divu;
-        IssmDouble   dvx[3],dvy[3],dvz[3];
         IssmDouble   *xyz_list = NULL;
         int numnodes = element->GetNumberOfNodes();
-
-        /*Initialize Element matrix and vectors*/
-        ElementVector* pe   = element->NewElementVector();
-        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -105 +100 @@
         element->GetVerticesCoordinates(&xyz_list);
 
+        /*Initialize Element matrix and vectors*/
+        ElementVector* pe    = element->NewElementVector();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dvx   = xNew<IssmDouble>(dim);
+        IssmDouble*    dvy   = xNew<IssmDouble>(dim);
+        IssmDouble*    dvz   = xNew<IssmDouble>(dim);
+
         Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
         Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-        Input* vz_input;
+        Input* vz_input = NULL;
         if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
 
@@ -122 +124 @@
                 }
 
+                divu=dvx[0]+dvy[1];
+                if (dim==3) divu=divu+dvz[2];
+
                 for(int i=0;i<numnodes;i++){
-                        divu=dvx[0]+dvy[1];
-                        if (dim==3) divu=divu+dvz[2];
-                        pe->values[i]+=divu*r*Jdet*gauss->weight*basis[i];
+                        pe->values[i] += - r * divu * Jdet * gauss->weight * basis[i];
                 }
         }
@@ -133 +136 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dvx);
+        xDelete<IssmDouble>(dvy);
+        xDelete<IssmDouble>(dvz);
         return pe;
 }/*}}}*/

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp

@@ -13 +13 @@
         Matrix<IssmDouble>*  Kff     = NULL;
         Matrix<IssmDouble>*  Kfs     = NULL;
-        Vector<IssmDouble>*  ug_old  = NULL;
         Vector<IssmDouble>*  ug      = NULL;
         Vector<IssmDouble>*  uf      = NULL;
@@ -23 +22 @@
         IssmDouble           eps_rel,r,theta; // 0<theta<.5   -> .15<theta<.45
         int                  configuration_type,max_nonlinear_iterations;
+        bool                 vel_converged      = false;
+        bool                 pressure_converged = false;
 
         /*Create analysis*/
@@ -40 +41 @@
         /*Convergence criterion*/
         int  count = 0;
-        GetSolutionFromInputsx(&ug,femmodel);
-        Vector<IssmDouble>* vx     = NULL;
-        Vector<IssmDouble>* vx_old = NULL;
-        GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
+        Vector<IssmDouble>* vel     = NULL;
+        Vector<IssmDouble>* vel_old = NULL;
+        GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+        GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
 
         while(true){
@@ -49 +50 @@
 
                 /*save pointer to old velocity*/
-                delete ug_old;ug_old=ug;
-                delete vx_old;vx_old=vx;
-                delete pug_old;pug_old=pug;
+                delete vel_old;vel_old=vel;
+                delete pug_old;pug_old=pug; 
+                pressure_converged=false; vel_converged=false;
 
                 /*Solve KU=F*/
                 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+                femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
                 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
                 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
@@ -63 +65 @@
 
                 /*Update solution*/
-                InputUpdateFromSolutionx(femmodel,ug); 
-                GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
+                InputUpdateFromSolutionx(femmodel,ug); delete ug;
+                GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
 
                 femmodel->SetCurrentConfiguration(UzawaPressureAnalysisEnum);
+                femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
                 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
                 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
@@ -72 +75 @@
                 Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
                 delete Kff; delete pf; delete df;
-                Mergesolutionfromftogx(&pug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+                Mergesolutionfromftogx(&pug, uf,ys,femmodel->nodes,femmodel->parameters); delete uf; delete ys;
 
                 /*Update solution*/
-                InputUpdateFromSolutionx(femmodel,pug); 
+                InputUpdateFromSolutionx(femmodel,pug); delete pug;
+                GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
 
                 /*Check for convergence*/
-                //Vector<IssmDouble>* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
-                //IssmDouble ndu=dug->Norm(NORM_TWO);   delete dug;
-                //IssmDouble nu =ug_old->Norm(NORM_TWO);
-                Vector<IssmDouble>* dvx=vx_old->Duplicate(); vx_old->Copy(dvx); dvx->AYPX(vx,-1.0);
-                IssmDouble ndu=dvx->Norm(NORM_TWO);   delete dvx;
-                IssmDouble nu =vx_old->Norm(NORM_TWO);
+                Vector<IssmDouble>* dvel=vel_old->Duplicate(); vel_old->Copy(dvel); dvel->AYPX(vel,-1.0);
+                IssmDouble ndu=dvel->Norm(NORM_TWO);   delete dvel;
+                IssmDouble nu =vel_old->Norm(NORM_TWO);
                 if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
                 if((ndu/nu)<eps_rel){
                         if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
-                        break;
+                        vel_converged=true;
                 }
                 else{ 
                         if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
+                        vel_converged=false;
+                }
+                Vector<IssmDouble>* dup=pug_old->Duplicate(); pug_old->Copy(dup); dup->AYPX(pug,-1.0);
+                IssmDouble ndp=dup->Norm(NORM_TWO);   delete dup;
+                IssmDouble np =pug_old->Norm(NORM_TWO);
+                if (xIsNan<IssmDouble>(ndp) || xIsNan<IssmDouble>(np)) _error_("convergence criterion is NaN!");
+                if((ndp/np)<eps_rel){
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(dp)/norm(p)" << ndp/np*100 << " < " << eps_rel*100 << " %\n");
+                        pressure_converged=true;
+                }
+                else{ 
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(dp/)/norm(p)" << ndp/np*100 << " > " << eps_rel*100 << " %\n");
+                        pressure_converged=false;
                 }
 
+                if(pressure_converged && vel_converged) break;
                 if(count>=max_nonlinear_iterations){
                         _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
@@ -101 +116 @@
         if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
 
-        delete ug;  
-        delete ug_old;  
         delete pug;  
         delete pug_old;  
-        delete vx;  
-        delete vx_old;  
+        delete vel;  
+        delete vel_old;  
         delete stressanalysis;
         delete pressureanalysis;