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

Stephane Popinet popinet at users.sf.net
Fri May 15 02:55:08 UTC 2009


The following commit has been merged in the upstream branch:
commit 3abf02a5cc6e92c2acc4c44288276165ae8c4317
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Thu Mar 13 10:25:11 2008 +1100

    Removed MRS acceleration
    
    While this accelerated convergence somewhat in some cases, it seemed to prevent
    convergence in other cases. This patch is not an exact "rollback" of the
    initial MRS implementation because it does not rollback other simple but
    important changes to the Poisson solver which seem to really improve robustness.
    
    darcs-hash:20080312232511-d4795-30d7551d50f18652854a7ad4b9a4fb087048583e.gz

diff --git a/src/event.c b/src/event.c
index c234066..a24d240 100644
--- a/src/event.c
+++ b/src/event.c
@@ -769,7 +769,7 @@ static void stream_from_vorticity (GfsDomain * domain,
   gfs_multilevel_params_init (&par);
   par.depth = gfs_domain_depth (domain);
   while (norm.infty > tolerance && maxit) {
-    gfs_poisson_cycle (domain, &par, stream, vorticity, dia, res, NULL, NULL);
+    gfs_poisson_cycle (domain, &par, stream, vorticity, dia, res);
     norm = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res);
     maxit--;
   }
diff --git a/src/ocean.c b/src/ocean.c
index af7c0a7..009122e 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -233,7 +233,7 @@ static void gfs_free_surface_pressure (GfsDomain * toplayer,
 	     par->residual.second, 
 	     par->residual.infty);
 #endif
-    gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1, NULL, NULL);
+    gfs_poisson_cycle (toplayer, par, p, fp.div, fp.dia, res1);
     par->residual = gfs_domain_norm_residual (toplayer, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
     par->niter++;
   }
diff --git a/src/poisson.c b/src/poisson.c
index fff8510..55b7d07 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -659,31 +659,6 @@ static void update_sv (FttCell * cell, MRSData * data)
   GFS_VALUE (cell, data->u) = GFS_VALUE (cell, data->v);
 }
 
-/* See Jun Zhang, "Multigrid acceleration techniques and applications
-   to the numerical solutions of partial differential equations", PhD
-   Thesis, George Washington University, section 4.2 */
-static void minimal_residual_smoothing (GfsDomain * domain, 
-					GfsVariable * v, GfsVariable * s, 
-					GfsVariable * u, GfsVariable * res)
-{
-  MRSData data;
-  data.s = s;
-  data.r = res;
-  data.u = u;
-  data.v = v;
-  data.srs = data.rs2 = 0.;
-  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			    (FttCellTraverseFunc) compute_beta, &data);
-  gfs_all_reduce (domain, data.srs, MPI_DOUBLE, MPI_SUM);
-  gfs_all_reduce (domain, data.rs2, MPI_DOUBLE, MPI_SUM);
-  if (data.rs2 > 0.) {
-    data.beta = data.srs/data.rs2;
-    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) update_sv, &data);
-    gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
-  }
-}
-
 /**
  * gfs_poisson_cycle:
  * @domain: the domain on which to solve the Poisson equation.
@@ -709,9 +684,7 @@ void gfs_poisson_cycle (GfsDomain * domain,
 			GfsVariable * u,
 			GfsVariable * rhs,
 			GfsVariable * dia,
-			GfsVariable * res,
-			GfsVariable * v,
-			GfsVariable * s)
+			GfsVariable * res)
 {
   guint n, l, nrelax, minlevel;
   GfsVariable * dp;
@@ -771,8 +744,6 @@ void gfs_poisson_cycle (GfsDomain * domain,
   gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, u);
   /* compute new residual on leaf cells */
   gfs_residual (domain, p->dimension, FTT_TRAVERSE_LEAFS, -1, u, rhs, dia, res);
-  if (v && s)
-    minimal_residual_smoothing (domain, v, s, u, res);
 
   gts_object_destroy (GTS_OBJECT (dp));
 }
diff --git a/src/poisson.h b/src/poisson.h
index 965fc5b..f26e5d9 100644
--- a/src/poisson.h
+++ b/src/poisson.h
@@ -73,10 +73,7 @@ void                  gfs_poisson_cycle              (GfsDomain * domain,
 						      GfsVariable * u,
 						      GfsVariable * rhs,
 						      GfsVariable * dia,
-						      GfsVariable * res,
-						      GfsVariable * v,
-						      GfsVariable * s);
-
+						      GfsVariable * res);
 void                  gfs_diffusion_coefficients     (GfsDomain * domain,
 						      GfsSourceDiffusion * d,
 						      gdouble dt,
diff --git a/src/simulation.c b/src/simulation.c
index 0ab1504..f70d13e 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1451,20 +1451,10 @@ static void copy_res (FttCell * cell, gpointer * data)
   GFS_VARIABLE (cell, res->i) = GFS_VARIABLE (cell, res1->i);
 }
 
-typedef struct {
-  GfsVariable * u, * r, * s, * v;
-} MRSData;
-
-static void init_mrs (FttCell * cell, MRSData * d)
-{
-  GFS_VALUE (cell, d->v) = GFS_VALUE (cell, d->u);
-  GFS_VALUE (cell, d->s) = GFS_VALUE (cell, d->r);
-}
-
 static void poisson_run (GfsSimulation * sim)
 {
   GfsDomain * domain = GFS_DOMAIN (sim);
-  GfsVariable * dia, * div, * res = NULL, * res1, * p, * s, * v;
+  GfsVariable * dia, * div, * res = NULL, * res1, * p;
   GfsMultilevelParams * par = &sim->approx_projection_params;
   GSList * i;
 
@@ -1484,18 +1474,11 @@ static void poisson_run (GfsSimulation * sim)
   gfs_poisson_coefficients (domain, NULL);
   res1 = gfs_temporary_variable (domain);
   dia = gfs_temporary_variable (domain);
-  s = gfs_temporary_variable (domain);
-  v = gfs_temporary_variable (domain);
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 			    (FttCellTraverseFunc) gfs_cell_reset, dia);
   /* compute residual */
   par->depth = gfs_domain_depth (domain);  
   gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
-  /* initialize Minimal Residual Smoothing */
-  MRSData mrs;
-  mrs.u = p; mrs.r = res1; mrs.v = v; mrs.s = s;
-  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-  			    (FttCellTraverseFunc) init_mrs, &mrs);
   /* solve for pressure */
   par->residual_before = par->residual = 
     gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
@@ -1517,7 +1500,7 @@ static void poisson_run (GfsSimulation * sim)
     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
 
     gfs_domain_timer_start (domain, "poisson_cycle");
-    gfs_poisson_cycle (domain, par, p, div, dia, res1, v, s);
+    gfs_poisson_cycle (domain, par, p, div, dia, res1);
     par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
     gfs_domain_timer_stop (domain, "poisson_cycle");
 
@@ -1541,8 +1524,6 @@ static void poisson_run (GfsSimulation * sim)
   gts_object_destroy (GTS_OBJECT (dia));
   gts_object_destroy (GTS_OBJECT (div));
   gts_object_destroy (GTS_OBJECT (res1));
-  gts_object_destroy (GTS_OBJECT (s));
-  gts_object_destroy (GTS_OBJECT (v));
 }
 
 static void poisson_class_init (GfsSimulationClass * klass)
diff --git a/src/timestep.c b/src/timestep.c
index 6c9d121..c60d883 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -272,16 +272,6 @@ void gfs_update_gradients (GfsDomain * domain,
   gfs_scale_gradients (domain, FTT_DIMENSION, g);
 }
 
-typedef struct {
-  GfsVariable * u, * r, * s, * v;
-} MRSData;
-
-static void init_mrs (FttCell * cell, MRSData * d)
-{
-  GFS_VALUE (cell, d->v) = GFS_VALUE (cell, d->u);
-  GFS_VALUE (cell, d->s) = GFS_VALUE (cell, d->r);
-}
-
 static void mac_projection (GfsDomain * domain,
 			    GfsMultilevelParams * par,
 			    GfsAdvectionParams * apar,
@@ -297,8 +287,6 @@ static void mac_projection (GfsDomain * domain,
   GfsVariable * div = gfs_temporary_variable (domain);
   GfsVariable * dia = gfs_temporary_variable (domain);
   GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
-  GfsVariable * s = gfs_temporary_variable (domain);
-  GfsVariable * v = gfs_temporary_variable (domain);
   /* Initialize face coefficients */
   gfs_poisson_coefficients (domain, alpha);
 
@@ -331,11 +319,6 @@ static void mac_projection (GfsDomain * domain,
   /* compute residual */
   par->depth = gfs_domain_depth (domain);
   gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
-  /* initialize Minimal Residual Smoothing */
-  MRSData mrs;
-  mrs.u = p; mrs.r = res1; mrs.v = v; mrs.s = s;
-  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-  			    (FttCellTraverseFunc) init_mrs, &mrs);
   /* solve for pressure */
   par->residual_before = par->residual = 
     gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
@@ -352,11 +335,7 @@ static void mac_projection (GfsDomain * domain,
 	     par->residual.second, 
 	     par->residual.infty);
 #endif
-#if 1
-    gfs_poisson_cycle (domain, par, p, div, dia, res1, v, s);
-#else
-    gfs_poisson_cycle (domain, par, p, div, dia, res1, NULL, NULL);
-#endif
+    gfs_poisson_cycle (domain, par, p, div, dia, res1);
     par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
     if (par->residual.infty == res_max_before) /* convergence has stopped!! */
       break;
@@ -372,8 +351,6 @@ static void mac_projection (GfsDomain * domain,
   gts_object_destroy (GTS_OBJECT (dia));
   if (!res)
     gts_object_destroy (GTS_OBJECT (res1));
-  gts_object_destroy (GTS_OBJECT (s));
-  gts_object_destroy (GTS_OBJECT (v));
 
   gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
   gfs_scale_gradients (domain, FTT_DIMENSION, g);

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list