[cig-commits] [commit] rajesh-petsc, rajesh-petsc-schur: Added calls to various PETSc versions of solve_Ahat_p_fhat functions (748ea48)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/citcoms

On branches: rajesh-petsc,rajesh-petsc-schur
Link       : https://github.com/geodynamics/citcoms/compare/464e1b32299b15819f93efd98d969cddb84dfe51...f97ae655a50bdbd6dac1923a3471ee4dae178fbd

>---------------------------------------------------------------

commit 748ea48a8ae9b2f523ae4bc1c45e2c54e7653cea
Author: Rajesh Kommu <rajesh.kommu at gmail.com>
Date:   Wed Aug 20 12:34:21 2014 -0700

    Added calls to various PETSc versions of solve_Ahat_p_fhat functions


>---------------------------------------------------------------

748ea48a8ae9b2f523ae4bc1c45e2c54e7653cea
 lib/Stokes_flow_Incomp.c | 59 ++++++++++++++++++++++++------------------------
 1 file changed, 30 insertions(+), 29 deletions(-)

diff --git a/lib/Stokes_flow_Incomp.c b/lib/Stokes_flow_Incomp.c
index 5038a1f..3a711a3 100644
--- a/lib/Stokes_flow_Incomp.c
+++ b/lib/Stokes_flow_Incomp.c
@@ -34,6 +34,7 @@
 #include <sys/types.h>
 #include "element_definitions.h"
 #include "global_defs.h"
+#include "petsc_citcoms.h"
 #include <stdlib.h>
 
 void myerror(struct All_variables *,char *);
@@ -93,27 +94,6 @@ void solve_constrained_flow_iterative(E)
 
     return;
 }
-/*
-PetscErrorCode solve_constrained_flow_iterative_petsc(struct All_variables *E)
-{
-  PetscErrorCode ierr;
-  int cycles;
-
-  cycles = E->control.p_iterations;
-  ierr = solve_Ahat_p_fhat_petsc(E, E->UVec, E->PVec, E->FVec, 
-      E->control.accuracy, &cycles);
-  CHKERRQ(ierr);
-
-  if(E->control.pseudo_free_surf)
-    v_from_vector_pseudo_surf_petsc(E);
-  else
-    v_from_vector_petsc(E);
-
-  p_to_nodes_petsc(E, E->PVec, E->NPVec, E->mesh.levmax);
-  
-  PetscFunctionReturn(0);
-}
-*/
 /* ========================================================================= */
 
 static double momentum_eqn_residual(struct All_variables *E,
@@ -201,6 +181,29 @@ static void solve_Ahat_p_fhat(struct All_variables *E,
                                double **V, double **P, double **F,
                                double imp, int *steps_max)
 {
+  if(E->control.use_petsc)
+  {
+    if(E->control.petsc_schur) // use Schur complement reduction
+    {
+      myerror(E, "Error: Schur complement reduction not implemented\n");
+    }
+    else                       // use the Uzawa algorithm
+    {
+      if(E->control.inv_gruneisen == 0)
+        solve_Ahat_p_fhat_CG_PETSc(E, V, P, F, imp, steps_max);
+      else
+      {
+        if(strcmp(E->control.uzawa, "cg") == 0)
+          solve_Ahat_p_fhat_iterCG(E, V, P, F, imp, steps_max);
+        else if(strcmp(E->control.uzawa, "bicg") == 0)
+          solve_Ahat_p_fhat_BiCG_PETSc(E, V, P, F, imp, steps_max);
+        else
+          myerror(E, "Error: unknown Uzawa iteration\n");
+      }
+    }
+  }
+  else                         // the original non-PETSc CitcomS code
+  {
     if(E->control.inv_gruneisen == 0)
         solve_Ahat_p_fhat_CG(E, V, P, F, imp, steps_max);
     else {
@@ -211,8 +214,7 @@ static void solve_Ahat_p_fhat(struct All_variables *E,
         else
             myerror(E, "Error: unknown Uzawa iteration\n");
     }
-
-    return;
+  }
 }
 
 
@@ -1204,15 +1206,10 @@ static void solve_Ahat_p_fhat_BiCG(struct All_variables *E,
         dpressure = sqrt(global_p_norm2(E, s0) / (1e-32 + E->monitor.pdotp));
 	
 
-	
-
         assemble_div_rho_u(E, V, t0, lev);
         E->monitor.incompressibility = sqrt(global_div_norm2(E, t0)
                                             / (1e-32 + E->monitor.vdotv));
 
-
-
-
         count++;
 
         if(E->control.print_convergence && E->parallel.me==0) {
@@ -1330,7 +1327,11 @@ static void solve_Ahat_p_fhat_iterCG(struct All_variables *E,
             for(i=1;i<=npno;i++) old_p[m][i] = P[m][i];
         }
 
-        solve_Ahat_p_fhat_CG(E, V, P, F, imp, &cycles);
+        if(E->control.use_petsc)
+          solve_Ahat_p_fhat_CG_PETSc(E, V, P, F, imp, &cycles);
+        else
+          solve_Ahat_p_fhat_CG(E, V, P, F, imp, &cycles);
+
 
         /* compute norm of div(rho*V) */
         assemble_div_rho_u(E, V, E->u1, lev);