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

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


The following commit has been merged in the upstream branch:
commit b910c6c919aa6f4be0cf243fe8be2bb38a61bdbb
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Fri Feb 22 08:35:36 2008 +1100

    Pressure term is included in RHS of viscosity solve
    
    This is an important change to the timestepping which should improve
    things significantly in particular when large source terms are included
    (e.g gravity + hydrostatic pressure).
    
    darcs-hash:20080221213536-d4795-6a7b40083fe3201219773bfa84f4d9ab315da60d.gz

diff --git a/configure.in b/configure.in
index 47262a1..9516ef2 100644
--- a/configure.in
+++ b/configure.in
@@ -12,8 +12,8 @@ dnl are available for $ac_help expansion (don't we all *love* autoconf?)
 # set GFS_BINARY_AGE and GFS_INTERFACE_AGE to 0.
 #
 GFS_MAJOR_VERSION=1
-GFS_MINOR_VERSION=1
-GFS_MICRO_VERSION=2
+GFS_MINOR_VERSION=2
+GFS_MICRO_VERSION=0
 GFS_INTERFACE_AGE=0
 GFS_BINARY_AGE=0
 GFS_VERSION=$GFS_MAJOR_VERSION.$GFS_MINOR_VERSION.$GFS_MICRO_VERSION
diff --git a/src/ocean.c b/src/ocean.c
index 7717024..312c684 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -111,24 +111,26 @@ static void gfs_correct_normal_velocities_weighted (GfsDomain * domain,
 						    gdouble dt,
 						    gboolean weighted)
 {
-  if (weighted)
-    gfs_correct_normal_velocities (domain, dimension, p, g, dt);
-  else {
-    gpointer data[3];
-    FttComponent c;
+  gpointer data[3];
+  FttComponent c;
     
-    g_return_if_fail (domain != NULL);
-    g_return_if_fail (p != NULL);
-    g_return_if_fail (g != NULL);
+  g_return_if_fail (domain != NULL);
+  g_return_if_fail (p != NULL);
+  g_return_if_fail (g != NULL);
     
-    for (c = 0; c < dimension; c++) {
-      g[c] = gfs_temporary_variable (domain);
-      gfs_variable_set_vector (g[c], c);
-    }
-    data[0] = g;
-    data[1] = &dimension;
-    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) reset_gradients, data);
+  for (c = 0; c < dimension; c++) {
+    g[c] = gfs_temporary_variable (domain);
+    gfs_variable_set_vector (g[c], c);
+  }
+  data[0] = g;
+  data[1] = &dimension;
+  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+			    (FttCellTraverseFunc) reset_gradients, data);
+  if (weighted) {
+    gfs_correct_normal_velocities (domain, dimension, p, g, dt);
+    gfs_scale_gradients (domain, dimension, g);
+  }
+  else {
     data[0] = p;
     data[1] = g;
     data[2] = &dt;
@@ -302,26 +304,17 @@ static void ocean_run (GfsSimulation * sim)
 
     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
 
-    gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., FALSE); 
-    /* just there so that the call below 
-       has sthg to free */
-    gfs_centered_velocity_advection_diffusion (domain, 2,
-					       &sim->advection_params,
-					       g,
-					       sim->physical_params.alpha);
-
     gfs_poisson_coefficients (domain, fH);
     gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., 
 					    sim->approx_projection_params.weighted);
-    if (gfs_has_source_coriolis (domain)) {
-      gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
-      gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
-      gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., 
-					      sim->approx_projection_params.weighted);
-      gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
-    }
-    else
-      gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+    gfs_centered_velocity_advection_diffusion (domain, 2,
+					       &sim->advection_params,
+					       g, g,
+					       sim->physical_params.alpha);
+    gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
+    gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
+    gts_object_destroy (GTS_OBJECT (g[0]));
+    gts_object_destroy (GTS_OBJECT (g[1]));
 
     sim->time.t = sim->tnext;
     sim->time.i++;
@@ -337,6 +330,8 @@ static void ocean_run (GfsSimulation * sim)
     gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2., 
 					    sim->approx_projection_params.weighted);
     gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+    gts_object_destroy (GTS_OBJECT (g[0]));
+    gts_object_destroy (GTS_OBJECT (g[1]));
     gfs_domain_timer_stop (domain, "free_surface_pressure");
 
     gfs_domain_cell_traverse (domain,
@@ -652,8 +647,8 @@ static void ocean_run (GfsSimulation * sim)
 
     gfs_domain_timer_start (domain, "correct_normal_velocities");
     gfs_poisson_coefficients (domain, NULL);
-    gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2., TRUE);
-
+    gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2.,
+					    sim->approx_projection_params.weighted);
     gfs_domain_cell_traverse_boundary (domain, FTT_BACK,
 				       FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 				       (FttCellTraverseFunc) compute_w, 
@@ -682,28 +677,19 @@ static void ocean_run (GfsSimulation * sim)
 
     gfs_centered_velocity_advection_diffusion (domain, 2,
 					       &sim->advection_params,
-					       g,
+					       g, g,
 					       sim->physical_params.alpha);
-
-    /* barotropic pressure and Coriolis terms */
-    gfs_poisson_coefficients (domain, NULL);
-    gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., 
-					    sim->approx_projection_params.weighted);
-    if (gfs_has_source_coriolis (domain)) {
-      gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
-      gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
-      gfs_correct_normal_velocities_weighted (domain, 2, p, g, 0., 
-					      sim->approx_projection_params.weighted);
-      gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
-    }
-    else
-      gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+    gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
+    gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt/2.);
 #else
     gfs_poisson_coefficients (domain, NULL);
     gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2.,
 					    sim->approx_projection_params.weighted);
     gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
 #endif
+    gts_object_destroy (GTS_OBJECT (g[0]));
+    gts_object_destroy (GTS_OBJECT (g[1]));
+
     sim->time.t = sim->tnext;
     sim->time.i++;
 
@@ -719,6 +705,9 @@ static void ocean_run (GfsSimulation * sim)
     gfs_correct_normal_velocities_weighted (domain, 2, p, g, sim->advection_params.dt/2.,
 					    sim->approx_projection_params.weighted);
     gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt/2.);
+    gts_object_destroy (GTS_OBJECT (g[0]));
+    gts_object_destroy (GTS_OBJECT (g[1]));
+    
     gfs_domain_timer_stop (domain, "free_surface_pressure");
 
     gfs_domain_cell_traverse (domain,
diff --git a/src/simulation.c b/src/simulation.c
index af03548..02462c3 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -332,7 +332,7 @@ static void advance_tracers (GfsDomain * domain, gdouble dt)
 
 static void simulation_run (GfsSimulation * sim)
 {
-  GfsVariable * p, * pmac, * res = NULL;
+  GfsVariable * p, * pmac, * res = NULL, * g[FTT_DIMENSION], * gmac[FTT_DIMENSION];
   GfsDomain * domain;
   GSList * i;
 
@@ -342,6 +342,13 @@ static void simulation_run (GfsSimulation * sim)
   g_assert (p);
   pmac = gfs_variable_from_name (domain->variables, "Pmac");
   g_assert (pmac);
+  FttComponent c;
+  for (c = 0; c < FTT_DIMENSION; c++) {
+    g[c] = gfs_temporary_variable (domain);
+    gmac[c] = gfs_temporary_variable (domain);
+    gfs_variable_set_vector (g[c], c);
+    gfs_variable_set_vector (gmac[c], c);
+  }
 
   gfs_simulation_refine (sim);
   gfs_simulation_init (sim);
@@ -358,13 +365,15 @@ static void simulation_run (GfsSimulation * sim)
     gfs_approximate_projection (domain,
 				&sim->approx_projection_params,
 				&sim->advection_params,
-				p, sim->physical_params.alpha, res);
+				p, sim->physical_params.alpha, res, g);
     gfs_simulation_set_timestep (sim);
     advance_tracers (domain, sim->advection_params.dt/2.);
   }
+  else
+    gfs_update_gradients (domain, p, sim->physical_params.alpha, g);
+
   while (sim->time.t < sim->time.end &&
 	 sim->time.i < sim->time.iend) {
-    GfsVariable * g[FTT_DIMENSION];
     gdouble tstart = gfs_clock_elapsed (domain->timer);
 
     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
@@ -375,7 +384,7 @@ static void simulation_run (GfsSimulation * sim)
     gfs_mac_projection (domain,
     			&sim->projection_params, 
     			&sim->advection_params,
-			p, sim->physical_params.alpha, g);
+			p, sim->physical_params.alpha, gmac);
     gfs_variables_swap (p, pmac);
 
     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
@@ -383,17 +392,11 @@ static void simulation_run (GfsSimulation * sim)
     gfs_centered_velocity_advection_diffusion (domain,
 					       FTT_DIMENSION,
 					       &sim->advection_params,
-					       g,
+					       gmac, sim->time.i > 0 ? g : gmac,
 					       sim->physical_params.alpha);
-
-    if (gfs_has_source_coriolis (domain)) {
-      gfs_poisson_coefficients (domain, sim->physical_params.alpha);
-      gfs_correct_normal_velocities (domain, 2, p, g, 0.);
-      gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);
-      gfs_source_coriolis_implicit (domain, sim->advection_params.dt);
-      gfs_correct_normal_velocities (domain, 2, p, g, 0.);
-      gfs_correct_centered_velocities (domain, 2, g, -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,
@@ -402,7 +405,7 @@ static void simulation_run (GfsSimulation * sim)
 
     gfs_approximate_projection (domain,
    				&sim->approx_projection_params, 
-    				&sim->advection_params, p, sim->physical_params.alpha, res);
+    				&sim->advection_params, p, sim->physical_params.alpha, res, g);
 
     sim->time.t = sim->tnext;
     sim->time.i++;
@@ -417,6 +420,11 @@ static void simulation_run (GfsSimulation * sim)
   }
   gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);  
   gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gts_object_destroy, NULL);
+
+  for (c = 0; c < FTT_DIMENSION; c++) {
+    gts_object_destroy (GTS_OBJECT (g[c]));
+    gts_object_destroy (GTS_OBJECT (gmac[c]));
+  }
 }
 
 static void gfs_simulation_class_init (GfsSimulationClass * klass)
@@ -1454,8 +1462,10 @@ static void poisson_run (GfsSimulation * sim)
   dia = gfs_temporary_variable (domain);
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 			    (FttCellTraverseFunc) gfs_cell_reset, dia);
-  par->depth = gfs_domain_depth (domain);
+  /* compute residual */
+  par->depth = gfs_domain_depth (domain);  
   gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
+  /* solve for pressure */
   par->residual_before = par->residual = 
     gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, 1., res1);
   par->niter = 0;
diff --git a/src/timestep.c b/src/timestep.c
index c75efe4..b4be6a0 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -25,7 +25,7 @@
 #include "solid.h"
 #include "tension.h"
 
-static void reset_gradients (FttCell * cell, gpointer * data)
+static void reset_cell_gradients (FttCell * cell, gpointer * data)
 {
   GfsVariable ** g = data[0];
   guint * dimension = data[1];    
@@ -35,6 +35,67 @@ static void reset_gradients (FttCell * cell, gpointer * data)
     GFS_VARIABLE (cell, g[c]->i) = 0.;
 }
 
+static void reset_gradients (GfsDomain * domain, guint dimension, GfsVariable ** g)
+{
+  gpointer data[2];
+  data[0] = g;
+  data[1] = &dimension;
+  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+			    (FttCellTraverseFunc) reset_cell_gradients, data);
+}
+
+static void scale_cell_gradients (FttCell * cell, gpointer * data)
+{
+  GfsVariable ** g = data[0];
+  guint * dimension = data[1];
+  FttComponent c;
+
+  if (GFS_IS_MIXED (cell)) {
+    GfsSolidVector * s = GFS_STATE (cell)->solid;
+
+    for (c = 0; c < *dimension; c++)
+      if (s->s[2*c] + s->s[2*c + 1] > 0.)
+	GFS_VARIABLE (cell, g[c]->i) /= s->s[2*c] + s->s[2*c + 1];
+      else
+	g_assert (GFS_VARIABLE (cell, g[c]->i) == 0.);
+  }
+  else {
+    FttCellNeighbors n;
+    
+    ftt_cell_neighbors (cell, &n);
+    for (c = 0; c < *dimension; c++) {
+      FttCell * c1 = n.c[2*c], * c2 = n.c[2*c + 1];
+      
+      if (c1 && c2 && !GFS_CELL_IS_GRADIENT_BOUNDARY (c1) && !GFS_CELL_IS_GRADIENT_BOUNDARY (c2))
+	GFS_VARIABLE (cell, g[c]->i) /= 2.;
+    }
+  }
+}
+
+/**
+ * gfs_scale_gradients:
+ * @domain: a #GfsDomain.
+ * @dimension: the number of dimensions.
+ * @g: the components of the gradient.
+ *
+ * Scales the gradient accumulated in @g (typically using
+ * gfs_correct_normal_velocities()).
+ */
+void gfs_scale_gradients (GfsDomain * domain, guint dimension, GfsVariable ** g)
+{
+  g_return_if_fail (domain != NULL);
+  g_return_if_fail (g != NULL);
+
+  gpointer data[2];
+  data[0] = g;
+  data[1] = &dimension;
+  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+			    (FttCellTraverseFunc) scale_cell_gradients, data);
+  FttComponent c;
+  for (c = 0; c < dimension; c++)
+    gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, g[c]);
+}
+
 static void correct_normal_velocity (FttCellFace * face,
 				     gpointer * data)
 {
@@ -70,34 +131,6 @@ static void correct_normal_velocity (FttCellFace * face,
     GFS_VARIABLE (face->neighbor, gv[c]->i) += dp*GFS_FACE_FRACTION_RIGHT (face);
 }
 
-static void scale_gradients (FttCell * cell, gpointer * data)
-{
-  GfsVariable ** g = data[0];
-  guint * dimension = data[1];
-  FttComponent c;
-
-  if (GFS_IS_MIXED (cell)) {
-    GfsSolidVector * s = GFS_STATE (cell)->solid;
-
-    for (c = 0; c < *dimension; c++)
-      if (s->s[2*c] + s->s[2*c + 1] > 0.)
-	GFS_VARIABLE (cell, g[c]->i) /= s->s[2*c] + s->s[2*c + 1];
-      else
-	g_assert (GFS_VARIABLE (cell, g[c]->i) == 0.);
-  }
-  else {
-    FttCellNeighbors n;
-    
-    ftt_cell_neighbors (cell, &n);
-    for (c = 0; c < *dimension; c++) {
-      FttCell * c1 = n.c[2*c], * c2 = n.c[2*c + 1];
-      
-      if (c1 && c2 && !GFS_CELL_IS_GRADIENT_BOUNDARY (c1) && !GFS_CELL_IS_GRADIENT_BOUNDARY (c2))
-	GFS_VARIABLE (cell, g[c]->i) /= 2.;
-    }
-  }
-}
-
 /**
  * gfs_correct_normal_velocities:
  * @domain: a #GfsDomain.
@@ -108,6 +141,9 @@ static void scale_gradients (FttCell * cell, gpointer * data)
  *
  * Corrects the normal velocity field of @domain using @p and and @dt.
  *
+ * Assumes that the Poisson weighting coefficients have already been
+ * computed using gfs_poisson_coefficients().
+ *
  * Also allocates the @g variables (if @g is not %NULL) and fills them
  * with the centered gradient of @p.
  */
@@ -118,35 +154,16 @@ void gfs_correct_normal_velocities (GfsDomain * domain,
 				    gdouble dt)
 {
   gpointer data[3];
-  FttComponent c;
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (p != NULL);
 
-  if (g) {
-    for (c = 0; c < dimension; c++) {
-      g[c] = gfs_temporary_variable (domain);
-      gfs_variable_set_vector (g[c], c);
-    }
-    data[0] = g;
-    data[1] = &dimension;
-    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) reset_gradients, data);
-  }
   data[0] = p;
   data[1] = g;
   data[2] = &dt;
   gfs_domain_face_traverse (domain, dimension == 2 ? FTT_XY : FTT_XYZ,
 			    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttFaceTraverseFunc) correct_normal_velocity, data);
-  if (g) {
-    data[0] = g;
-    data[1] = &dimension;
-    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) scale_gradients, data);
-    for (c = 0; c < dimension; c++)
-      gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, g[c]);
-  }
 }
 
 static void scale_divergence (FttCell * cell, gpointer * data)
@@ -173,8 +190,6 @@ static void surface_tension (GfsDomain * domain,
 	gfs_correct_normal_velocities (domain, FTT_DIMENSION,
 				       GFS_SOURCE_TENSION_GENERIC (s)->c,
 				       g, dt);
-	if (g)
-	  gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, dt);
       }
       i = i->next;
     }
@@ -182,66 +197,59 @@ static void surface_tension (GfsDomain * domain,
 }
 
 /**
- * gfs_mac_projection:
+ * gfs_update_gradients:
  * @domain: a #GfsDomain.
- * @par: the projection control parameters.
- * @apar: the advection parameters.
  * @p: the pressure.
  * @alpha: the Poisson equation gradient weight.
  * @g: where to store the pressure gradient.
  *
- * Corrects the face-centered velocity field (MAC field) on the leaf
- * level of @domain using an exact (MAC) projection. The resulting
- * face-centered velocity field is (almost) exactly divergence
- * free. The (potential) pressure field is also obtained as a
- * by-product as well as its gradient at the center of the leaf cells
- * of the domain. The gradient is stored in newly-allocated @g[]
- * variables and is obtained by simple averaging from the face values
- * to the center. The newly-allocated @g[] variables should be freed
- * when not needed anymore.
- *
- * The @residual field of the @par projection parameters is set to the
- * norm of the residual after the projection. The @niter field of the
- * @par projection parameters is set to the number of iterations
- * performed to solve the Poisson equation. The other projection
- * parameters are not modified.
+ * Updates the gradients in @g using @p and @alpha.
  */
-void gfs_mac_projection (GfsDomain * domain,
-			 GfsMultilevelParams * par,
-			 GfsAdvectionParams * apar,
-			 GfsVariable * p,
-			 GfsFunction * alpha,
-			 GfsVariable ** g)
+void gfs_update_gradients (GfsDomain * domain, 
+			   GfsVariable * p,  
+			   GfsFunction * alpha, 
+			   GfsVariable ** g)
 {
-  gdouble dt;
-  gpointer data[2];
-  GfsVariable * div, * dia, * res;
-
   g_return_if_fail (domain != NULL);
-  g_return_if_fail (par != NULL);
-  g_return_if_fail (apar != NULL);
   g_return_if_fail (p != NULL);
   g_return_if_fail (g != NULL);
 
-  gfs_domain_timer_start (domain, "mac_projection");
-
-  div = gfs_temporary_variable (domain);
-  dia = gfs_temporary_variable (domain);
-  res = gfs_temporary_variable (domain);
-  
-  apar->v = gfs_variable_from_name (domain->variables, "U");
-  dt = apar->dt;
-  apar->dt /= 2.;
+  /* Add surface tension */
+  reset_gradients (domain, FTT_DIMENSION, g);
+  surface_tension (domain, gfs_domain_velocity (domain)[0], 0., alpha, g);
+  /* Initialize face coefficients */
+  gfs_poisson_coefficients (domain, alpha);
+  /* Add pressure gradient */
+  gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, 0.);
+  gfs_scale_gradients (domain, FTT_DIMENSION, g);
+}
 
+static void mac_projection (GfsDomain * domain,
+			    GfsMultilevelParams * par,
+			    GfsAdvectionParams * apar,
+			    GfsVariable * p,
+			    GfsFunction * alpha,
+			    GfsVariable * res,
+			    GfsVariable ** g)
+{
   /* Add surface tension */
-  surface_tension (domain, apar->v, apar->dt, alpha, NULL);
+  reset_gradients (domain, FTT_DIMENSION, g);
+  surface_tension (domain, gfs_domain_velocity (domain)[0], apar->dt, alpha, g);
+
+  GfsVariable * div = gfs_temporary_variable (domain);
+  GfsVariable * dia = gfs_temporary_variable (domain);
+  GfsVariable * res1 = res ? res : gfs_temporary_variable (domain);
   /* Initialize face coefficients */
   gfs_poisson_coefficients (domain, alpha);
+
+  /* Initialize diagonal coefficient */
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 			    (FttCellTraverseFunc) gfs_cell_reset, dia);
+
   /* compute MAC divergence */
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttCellTraverseFunc) gfs_normal_divergence, div);
+  gpointer data[2];
   data[0] = div;
   data[1] = &apar->dt;
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -259,25 +267,84 @@ void gfs_mac_projection (GfsDomain * domain,
 	     norm.first, norm.second, norm.infty);
   }
 #endif
-
-  /* solve for pressure */
+  
+  /* compute residual */
   par->depth = gfs_domain_depth (domain);
-  gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res);
+  gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
+  /* solve for pressure */
   par->residual_before = par->residual = 
-    gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
+    gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
   par->niter = 0;
   while (par->niter < par->nitermin ||
 	 (par->residual.infty > par->tolerance && par->niter < par->nitermax)) {
-    gfs_poisson_cycle (domain, par, p, div, dia, res);
-    par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res);
+#if 0
+    fprintf (stderr, "%d bias: %g first: %g second: %g infty: %g\n",
+	     par->niter, 
+	     par->residual.bias, 
+	     par->residual.first, 
+	     par->residual.second, 
+	     par->residual.infty);
+#endif
+    gfs_poisson_cycle (domain, par, p, div, dia, res1);
+    par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
     par->niter++;
   }
 
   gts_object_destroy (GTS_OBJECT (div));
   gts_object_destroy (GTS_OBJECT (dia));
-  gts_object_destroy (GTS_OBJECT (res));
+  if (!res)
+    gts_object_destroy (GTS_OBJECT (res1));
 
   gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
+  gfs_scale_gradients (domain, FTT_DIMENSION, g);
+}
+
+/**
+ * gfs_mac_projection:
+ * @domain: a #GfsDomain.
+ * @par: the projection control parameters.
+ * @apar: the advection parameters.
+ * @p: the pressure.
+ * @alpha: the Poisson equation gradient weight.
+ * @g: where to store the pressure gradient.
+ *
+ * Corrects the face-centered velocity field (MAC field) on the leaf
+ * level of @domain using an exact (MAC) projection. The resulting
+ * face-centered velocity field is (almost) exactly divergence
+ * free. The (potential) pressure field is also obtained as a
+ * by-product as well as its gradient at the center of the leaf cells
+ * of the domain. The gradient is stored in newly-allocated @g[]
+ * variables and is obtained by simple averaging from the face values
+ * to the center. The newly-allocated @g[] variables should be freed
+ * when not needed anymore.
+ *
+ * The @residual field of the @par projection parameters is set to the
+ * norm of the residual after the projection. The @niter field of the
+ * @par projection parameters is set to the number of iterations
+ * performed to solve the Poisson equation. The other projection
+ * parameters are not modified.
+ */
+void gfs_mac_projection (GfsDomain * domain,
+			 GfsMultilevelParams * par,
+			 GfsAdvectionParams * apar,
+			 GfsVariable * p,
+			 GfsFunction * alpha,
+			 GfsVariable ** g)
+{
+  gdouble dt;
+
+  g_return_if_fail (domain != NULL);
+  g_return_if_fail (par != NULL);
+  g_return_if_fail (apar != NULL);
+  g_return_if_fail (p != NULL);
+  g_return_if_fail (g != NULL);
+
+  gfs_domain_timer_start (domain, "mac_projection");
+
+  dt = apar->dt;
+  apar->dt /= 2.;
+
+  mac_projection (domain, par, apar, p, alpha, NULL, g);
 
 #if 0
   {
@@ -337,11 +404,8 @@ void gfs_correct_centered_velocities (GfsDomain * domain,
   data[3] = &dimension;
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttCellTraverseFunc) correct, data);
-  for (c = 0; c < dimension; c++) {
+  for (c = 0; c < dimension; c++)
     gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, v[c]);
-    gts_object_destroy (GTS_OBJECT (g[c]));
-    g[c] = NULL;
-  }
 }
 
 /**
@@ -375,15 +439,14 @@ void gfs_approximate_projection (GfsDomain * domain,
 				 GfsAdvectionParams * apar,
 				 GfsVariable * p,
 				 GfsFunction * alpha,
-				 GfsVariable * res)
+				 GfsVariable * res,
+				 GfsVariable ** g)
 {
-  gpointer data[2];
-  GfsVariable * dia, * div, * g[FTT_DIMENSION], * res1;
-
   g_return_if_fail (domain != NULL);
   g_return_if_fail (par != NULL);
   g_return_if_fail (apar != NULL);
   g_return_if_fail (p != NULL);
+  g_return_if_fail (g != NULL);
 
   gfs_domain_timer_start (domain, "approximate_projection");
   
@@ -395,67 +458,9 @@ void gfs_approximate_projection (GfsDomain * domain,
 			    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttFaceTraverseFunc) gfs_face_interpolated_normal_velocity, 
 			    gfs_domain_velocity (domain));
-  /* Add surface tension */
-  surface_tension (domain, gfs_domain_velocity (domain)[0], apar->dt, alpha, g);
-
-  div = gfs_temporary_variable (domain);
-  dia = gfs_temporary_variable (domain);
-  res1 = res ? res : gfs_temporary_variable (domain);
-  /* Initialize face coefficients */
-  gfs_poisson_coefficients (domain, alpha);
-
-  /* Initialize diagonal coefficient */
-  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
-			    (FttCellTraverseFunc) gfs_cell_reset, dia);
-
-  /* compute MAC divergence */
-  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			    (FttCellTraverseFunc) gfs_normal_divergence, div);
-  data[0] = div;
-  data[1] = &apar->dt;
-  gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-  			    (FttCellTraverseFunc) scale_divergence, data);
-
-#if 0
-  {
-    FILE * fp = fopen ("/tmp/macapprox", "wt");
-    GfsNorm norm;
-
-    gfs_write_mac_velocity (domain, 0.9, FTT_TRAVERSE_LEAFS, -1, NULL, fp);
-    fclose (fp);
-    norm = gfs_domain_norm_variable (domain, div, FTT_TRAVERSE_LEAFS, -1);
-    fprintf (stderr, "mac div before: %g %g %g\n",
-	     norm.first, norm.second, norm.infty);
-  }
-#endif
   
-  /* solve for pressure */
-  par->depth = gfs_domain_depth (domain);
-  gfs_residual (domain, par->dimension, FTT_TRAVERSE_LEAFS, -1, p, div, dia, res1);
-  par->residual_before = par->residual = 
-    gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
-  par->niter = 0;
-  while (par->niter < par->nitermin ||
-	 (par->residual.infty > par->tolerance && par->niter < par->nitermax)) {
-#if 0
-    fprintf (stderr, "%d bias: %g first: %g second: %g infty: %g\n",
-	     par->niter, 
-	     par->residual.bias, 
-	     par->residual.first, 
-	     par->residual.second, 
-	     par->residual.infty);
-#endif
-    gfs_poisson_cycle (domain, par, p, div, dia, res1);
-    par->residual = gfs_domain_norm_residual (domain, FTT_TRAVERSE_LEAFS, -1, apar->dt, res1);
-    par->niter++;
-  }
+  mac_projection (domain, par, apar, p, alpha, res, g);
 
-  gts_object_destroy (GTS_OBJECT (div));
-  gts_object_destroy (GTS_OBJECT (dia));
-  if (!res)
-    gts_object_destroy (GTS_OBJECT (res1));
-
-  gfs_correct_normal_velocities (domain, FTT_DIMENSION, p, g, apar->dt);
   gfs_correct_centered_velocities (domain, FTT_DIMENSION, g, apar->dt);
 
   gfs_domain_timer_stop (domain, "approximate_projection");
@@ -582,16 +587,22 @@ static GfsSourceDiffusion * source_diffusion (GfsVariable * v)
   return NULL;
 }
 
+static void add_pressure_gradient (FttCell * cell, GfsAdvectionParams * par)
+{
+  GFS_VALUE (cell, par->fv) -= GFS_VALUE (cell, par->g[par->v->component])*par->dt;
+}
+
 static void variable_sources (GfsDomain * domain,
 			      GfsAdvectionParams * par,
 			      GfsVariable * sv,
+			      GfsVariable ** gmac,
 			      GfsVariable ** g)
 {
   if (par->scheme == GFS_GODUNOV) {
     GfsVariable * v = par->v;
 
     par->u = gfs_domain_velocity (domain);
-    par->g = g;
+    par->g = gmac;
     par->fv = gfs_temporary_variable (domain);
     par->upwinding = GFS_FACE_UPWINDING;
     gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -609,6 +620,14 @@ static void variable_sources (GfsDomain * domain,
     gts_object_destroy (GTS_OBJECT (par->fv));
     par->fv = NULL;
   }
+  if (g) {
+    par->fv = sv;
+    par->g = g;
+    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+			      (FttCellTraverseFunc) add_pressure_gradient, par);
+    par->g = NULL;
+    par->fv = NULL;
+  }
   /* fixme: time should be set to t + dt/2 here for evaluation of
      source terms in the call below */
   gfs_domain_variable_centered_sources (domain, par->v, sv, par->dt);
@@ -640,6 +659,7 @@ static void variable_diffusion (GfsDomain * domain,
  * @domain: a #GfsDomain.
  * @dimension: the number of dimensions (2 or 3).
  * @apar: the advection parameters.
+ * @gmac: the MAC pressure gradient.
  * @g: the pressure gradient.
  * @alpha: the inverse of density or %NULL.
  *
@@ -659,6 +679,7 @@ static void variable_diffusion (GfsDomain * domain,
 void gfs_centered_velocity_advection_diffusion (GfsDomain * domain,
 						guint dimension,
 						GfsAdvectionParams * apar,
+						GfsVariable ** gmac,
 						GfsVariable ** g,
 						GfsFunction * alpha)
 {
@@ -667,7 +688,7 @@ void gfs_centered_velocity_advection_diffusion (GfsDomain * domain,
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (apar != NULL);
-  g_return_if_fail (g != NULL);
+  g_return_if_fail (gmac != NULL);
 
   gfs_domain_timer_start (domain, "centered_velocity_advection_diffusion");
 
@@ -683,20 +704,15 @@ void gfs_centered_velocity_advection_diffusion (GfsDomain * domain,
       rhs = gfs_temporary_variable (domain);
       gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 				(FttCellTraverseFunc) gfs_cell_reset, rhs);
-      variable_sources (domain, apar, rhs, g);
-      gts_object_destroy (GTS_OBJECT (g[c]));
-      g[c] = NULL;
+      variable_sources (domain, apar, rhs, gmac, g);
       variable_diffusion (domain, d, apar, rhs, alpha);
       gts_object_destroy (GTS_OBJECT (rhs));
     }
     else {
-      variable_sources (domain, apar, apar->v, g);
-      gts_object_destroy (GTS_OBJECT (g[c]));
-      g[c] = NULL;
+      variable_sources (domain, apar, apar->v, gmac, g);
       gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, apar->v);
     }
   }
-
   gfs_domain_timer_stop (domain, "centered_velocity_advection_diffusion");
 }
 
@@ -724,12 +740,12 @@ void gfs_tracer_advection_diffusion (GfsDomain * domain,
     rhs = gfs_temporary_variable (domain);
     gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			      (FttCellTraverseFunc) gfs_cell_reset, rhs);
-    variable_sources (domain, par, rhs, NULL);
+    variable_sources (domain, par, rhs, NULL, NULL);
     variable_diffusion (domain, d, par, rhs, NULL);
     gts_object_destroy (GTS_OBJECT (rhs));
   }
   else {
-    variable_sources (domain, par, par->v, NULL);
+    variable_sources (domain, par, par->v, NULL, NULL);
     gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, par->v);
   }
 
diff --git a/src/timestep.h b/src/timestep.h
index 5985bc9..db11a99 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -33,6 +33,13 @@ void          gfs_correct_normal_velocities   (GfsDomain * domain,
 					       GfsVariable * p,
 					       GfsVariable ** g,
 					       gdouble dt);
+void          gfs_scale_gradients             (GfsDomain * domain, 
+					       guint dimension, 
+					       GfsVariable ** g);
+void          gfs_update_gradients            (GfsDomain * domain, 
+					       GfsVariable * p,  
+					       GfsFunction * alpha, 
+					       GfsVariable ** g);
 void          gfs_mac_projection              (GfsDomain * domain,
 					       GfsMultilevelParams * par,
 					       GfsAdvectionParams * apar,
@@ -48,7 +55,8 @@ void          gfs_approximate_projection      (GfsDomain * domain,
 					       GfsAdvectionParams * apar,
 					       GfsVariable * p,
 					       GfsFunction * alpha,
-					       GfsVariable * res);
+					       GfsVariable * res,
+					       GfsVariable ** g);
 void          gfs_predicted_face_velocities   (GfsDomain * domain,
 					       guint d,
 					       GfsAdvectionParams * par);
@@ -61,6 +69,7 @@ void          gfs_diffusion                   (GfsDomain * domain,
 void          gfs_centered_velocity_advection_diffusion (GfsDomain * domain,
 							 guint dimension,
 							 GfsAdvectionParams * apar,
+							 GfsVariable ** gmac,
 							 GfsVariable ** g,
 							 GfsFunction * alpha);
 void          gfs_tracer_advection_diffusion  (GfsDomain * domain,

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list