Programming the Advection Scheme

From Gerris

The core of a fluid simulation

Let us assume that as a computational physicists you are most interested in how the Navier-Stokes equations are read. You are familiar with Stephane's J. Comput. Phys. article and you want to see how it is implemented in code. The core of the simulation in simulation.c specifically the following lines. We start with some bookkeeping:

static void simulation_run (GfsSimulation * sim)
{
  GfsVariable * p, * pmac, * res = NULL;
  GfsDomain * domain;
  GSList * i;
 
  domain = GFS_DOMAIN (sim);
 
  p = gfs_variable_from_name (domain->variables, "P");
  g_assert (p);
  pmac = gfs_variable_from_name (domain->variables, "Pmac");
  g_assert (pmac);
 
  gfs_simulation_refine (sim);
  gfs_simulation_init (sim);
 
  i = domain->variables;
  while (i) {
    if (GFS_IS_VARIABLE_RESIDUAL (i->data))
      res = i->data;
    i = i->next;
  }

gfs_simulation_set_timestep (sim);
  if (sim->time.i == 0) {
    gfs_approximate_projection (domain,
				&sim->approx_projection_params,
				&sim->advection_params,
				p, sim->physical_params.alpha, res);
    advance_tracers (domain, sim->advection_params.dt/2.);
  }
  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);
 
    gfs_simulation_set_timestep (sim);
 
    gfs_predicted_face_velocities (domain, FTT_DIMENSION, &sim->advection_params);
 
    gfs_variables_swap (p, pmac);
    gfs_mac_projection (domain,
    			&sim->projection_params, 
    			&sim->advection_params,
			p, sim->physical_params.alpha, g);
    gfs_variables_swap (p, pmac);
 
    gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);
 
    gfs_centered_velocity_advection_diffusion (domain,
					       FTT_DIMENSION,
					       &sim->advection_params,
					       g,
					       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_domain_cell_traverse (domain,
			      FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
			      (FttCellTraverseFunc) gfs_cell_coarse_init, domain);
    gfs_simulation_adapt (sim);
 
    gfs_approximate_projection (domain,
   				&sim->approx_projection_params, 
    				&sim->advection_params, p, sim->physical_params.alpha, res);
 
    advance_tracers (domain, sim->advection_params.dt);
 
    sim->time.t = sim->tnext;
    sim->time.i++;
 
    gts_range_add_value (&domain->timestep, gfs_clock_elapsed (domain->timer) - tstart);
    gts_range_update (&domain->timestep);
    gts_range_add_value (&domain->size, gfs_domain_size (domain, FTT_TRAVERSE_LEAFS, -1));
    gts_range_update (&domain->size);
  }
  gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);  
  gts_container_foreach (GTS_CONTAINER (sim->events),
			 (GtsFunc) gts_object_destroy, NULL);
}