[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 96204e7d10e8f81937972d250752c7a02cdbbc8b
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Wed Feb 27 08:32:11 2008 +1100

    Pressure-correction in the Crank-Nicholson scheme is turned off by default
    
    This was degrading the results for a number of test cases (capillary waves in
    particular). The reasons why are still unclear.
    
    Note that this does not mean that this version of the code is the same as the
    version prior to the "Pressure term is included in RHS of viscosity solve"
    patch, as this version also includes an important bug fix for the advection
    terms when surface tension is present.
    
    The pressure-correction Crank-Nicholson scheme can be turned off explicitly
    when required (e.g. "hydrostatic" test case).
    
    darcs-hash:20080226213211-d4795-81bcff59de78bc33c77dc942cc8305c601942add.gz

diff --git a/src/advection.c b/src/advection.c
index cb3d0b0..25df4e5 100644
--- a/src/advection.c
+++ b/src/advection.c
@@ -921,6 +921,8 @@ void gfs_advection_params_write (GfsAdvectionParams * par, FILE * fp)
 	   par->flux == gfs_face_velocity_convective_flux ?
 	   "gfs_face_velocity_convective_flux" : "NULL",
 	   par->average);
+  if (par->gc)
+    fputs ("  gc       = 1\n", fp);
   switch (par->scheme) {
   case GFS_GODUNOV: fputs ("  scheme   = godunov\n", fp); break;
   case GFS_NONE:    fputs ("  scheme   = none\n", fp); break;
@@ -942,6 +944,7 @@ void gfs_advection_params_init (GfsAdvectionParams * par)
   par->use_centered_velocity = TRUE;
   par->scheme = GFS_GODUNOV;
   par->average = FALSE;
+  par->gc = FALSE;
 }
 
 void gfs_advection_params_read (GfsAdvectionParams * par, GtsFile * fp)
@@ -952,6 +955,7 @@ void gfs_advection_params_read (GfsAdvectionParams * par, GtsFile * fp)
     {GTS_STRING, "flux",     TRUE},
     {GTS_STRING, "scheme",   TRUE},
     {GTS_INT,    "average",  TRUE},
+    {GTS_INT,    "gc",       TRUE},
     {GTS_NONE}
   };
   gchar * gradient = NULL, * flux = NULL, * scheme = NULL;
@@ -964,6 +968,7 @@ void gfs_advection_params_read (GfsAdvectionParams * par, GtsFile * fp)
   var[2].data = &flux;
   var[3].data = &scheme;
   var[4].data = &par->average;
+  var[5].data = &par->gc;
 
   gts_file_assign_variables (fp, var);
 
diff --git a/src/advection.h b/src/advection.h
index 17e6f13..32cd468 100644
--- a/src/advection.h
+++ b/src/advection.h
@@ -50,7 +50,7 @@ struct _GfsAdvectionParams {
   GfsUpwinding upwinding;
   GfsFaceAdvectionFluxFunc flux;
   GfsAdvectionScheme scheme;
-  gboolean average;
+  gboolean average, gc;
 };
 
 void         gfs_advection_params_init        (GfsAdvectionParams * par);
diff --git a/src/simulation.c b/src/simulation.c
index 02462c3..5c8c6c4 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -333,6 +333,7 @@ static void advance_tracers (GfsDomain * domain, gdouble dt)
 static void simulation_run (GfsSimulation * sim)
 {
   GfsVariable * p, * pmac, * res = NULL, * g[FTT_DIMENSION], * gmac[FTT_DIMENSION];
+  GfsVariable ** gc = sim->advection_params.gc ? g : NULL;
   GfsDomain * domain;
   GSList * i;
 
@@ -392,11 +393,21 @@ static void simulation_run (GfsSimulation * sim)
     gfs_centered_velocity_advection_diffusion (domain,
 					       FTT_DIMENSION,
 					       &sim->advection_params,
-					       gmac, sim->time.i > 0 ? g : gmac,
+					       gmac,
+					       sim->time.i > 0 || !gc ? gc : gmac,
 					       sim->physical_params.alpha);
-    gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
-    gfs_correct_centered_velocities (domain, FTT_DIMENSION, sim->time.i > 0 ? g : gmac, 
-				     -sim->advection_params.dt);
+    if (gc) {
+      gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
+      gfs_correct_centered_velocities (domain, FTT_DIMENSION, sim->time.i > 0 ? gc : gmac, 
+				       -sim->advection_params.dt);
+    }
+    else if (gfs_has_source_coriolis (domain)) {
+      gfs_correct_centered_velocities (domain, FTT_DIMENSION, sim->time.i > 0 ? g : gmac,
+				       sim->advection_params.dt);
+      gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
+      gfs_correct_centered_velocities (domain, FTT_DIMENSION, sim->time.i > 0 ? g : gmac, 
+				       -sim->advection_params.dt);
+    }
 
     gfs_domain_cell_traverse (domain,
 			      FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
diff --git a/test/hydrostatic/hydrostatic.gfs b/test/hydrostatic/hydrostatic.gfs
index 32e9c82..8090ea9 100644
--- a/test/hydrostatic/hydrostatic.gfs
+++ b/test/hydrostatic/hydrostatic.gfs
@@ -17,6 +17,11 @@
     Time { iend = 10 }
     ApproxProjectionParams { tolerance = 1e-12 }
     ProjectionParams { tolerance = 1e-12 }
+
+    # The pressure correction needs to be included in the
+    # Crank-Nicholson scheme in order to balance the body force
+    AdvectionParams { gc = 1 }
+
     OutputScalarNorm { istep = 1 } v { v = V }
     EventScript { start = end } {
         if awk '{if ($9 > 1e-12) exit (1);}' < v ; then

-- 
Gerris Flow Solver