7.1
general documentation
Groundwater flow module using CDO schemes

The groundwater flow module (GWF) relies on CDO vertex-based or CDO vertex+cell-based schemes.

To set-up a GWF computation, one has to update the cs_user_parameters.c file and edit the function cs_user_model at least in simple cases. In more complex cases, editing cs_user_finalize_setup should be necessary.

Activate the GWF module

The first step is to activate the CDO module in the function cs_user_model (please refer to Activation of CDO/HHO schemes).

Then, one has to activate the groundwater flow module (GWF) in the function cs_user_model. The function to call is cs_gwf_activate.

There are three parameters:

  1. The type of permeability to consider. Two choices: CS_PROPERTY_ISO or CS_PROPERTY_ANISO
  2. The main model to consider (i.e. which equations to solve). Please refer to Main models available in the GWF
  3. Optional flags to specify a physical phenomena to take into account or to specify a numerical treatment to apply.

Here are listed the available option flags:

Main models available in the GWF

There are currently two models :

First example: Activate the GWF model with an isotropic permeability and a fully saturated single-phase flow model.

{
/* For the groundwater flow module:
cs_gwf_activate(permeability_type,
model_type,
option_flag);
If option_flag = 0, then there is no option to set.
*/
cs_flag_t option_flag = 0;
cs_gwf_activate(CS_PROPERTY_ISO,
CS_GWF_MODEL_SATURATED_SINGLE_PHASE,
option_flag);
}

Second example: Activate the GWF model with an anisotropic permeability and an unsaturated single-phase flow model. Moreover, one takes into account the gravity effect.

{
/* Take into account the gravity effect */
cs_flag_t option_flag = CS_GWF_GRAVITATION;
/* In this case, the gravity vector has to be defined (either using the GUI
or in cs_user_parameters() function */
cs_gwf_activate(CS_PROPERTY_ANISO,
CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE,
option_flag);
}

Soils

The second step is to add at least one new soil. The add of soil(s) should be done before adding tracers. The function which adds a soil is cs_gwf_add_soil to be located inside the function cs_user_model. The last parameter of this function is the bulk density. This parameter is only useful if there is at least one tracer defined.

Predefined soil models

There are two predefined soil models

and one user-defined soil model

Remark: If a CS_GWF_MODEL_SATURATED_SINGLE_PHASE is set at the activation step, then one expects that all soil models are defined by the type CS_GWF_SOIL_SATURATED

Case of a saturated model

If the soil is saturated then this is a simple definition which can be located inside the function cs_user_model.

Example for two saturated soils defined by an anisotropic (saturated) permeability

{
/* Example 1: two "saturated" and anisotropic soils */
/* saturated anisotropic permeability */
const double aniso_val1[3][3] = {{1e-3, 0, 0},
{ 0, 1, 0},
{ 0, 0, 1e-1}};
const double theta_s = 1;
const cs_real_t bulk_density = 1.0; /* useless if no tracer is
considered */
/* Two (volume) zones have be defined called "soil1" and "soil2". */
cs_gwf_soil_t *s1 = cs_gwf_add_soil("soil1",
bulk_density,
theta_s,
CS_GWF_SOIL_SATURATED);
cs_gwf_soil_set_aniso_saturated(s1, aniso_val1);
/* For "simple" soil definitions, definition can be made here. */
cs_gwf_soil_t *s2 = cs_gwf_add_soil("soil2",
bulk_density,
theta_s,
CS_GWF_SOIL_SATURATED);
/* saturated anisotropic permeability */
const double aniso_val2[3][3] = {{1e-5, 0, 0},
{ 0, 1, 0},
{ 0, 0, 1e-2}};
cs_gwf_soil_set_aniso_saturated(s2, aniso_val2);
}

Use cs_gwf_set_iso_saturated_soil to set a saturated soil defined by an isotropic permeability instead of an anisotropic permeability. In this case, the second parameter is simply a ``double''.

Case of a Van Genuchten-Mualen model

Soils which behave according to a Van Genuchten-Mualen model can be specified using cs_gwf_soil_set_aniso_genuchten or cs_gwf_soil_set_iso_genuchten according to the type of permeability which has been defined.

Example of a soil relying on a Van Genuchten-Mualen and considering a isotropic permeability

{
/* Example 2: Add a new user-defined soil for all the cells */
const cs_real_t theta_s = 0.9; /* max. liquid saturation */
const cs_real_t bulk_density = 1800.0; /* useless if no tracer is
considered */
cs_gwf_soil_t *s = cs_gwf_add_soil("cells", /* volume zone name */
bulk_density,
theta_s,
CS_GWF_SOIL_GENUCHTEN); /* soil model */
cs_gwf_soil_set_iso_genuchten(s,
3e-1, /* saturated permeability */
0.078, /* residual moisture */
0.036, /* scaling parameter */
1.56, /* (n) shape parameter */
0.5); /* (L) tortuosity */
}

User-defined soils

If the predefined models are not sufficient, it is possible to add a user-defined soil. In this case, the add of the new soil is made as follows

  1. Add a user-defined model and set constant members of the soil structure and then call cs_gwf_soil_set_user (inside the function cs_user_model)
  2. Set the remaing model parameters (inside the function cs_user_finalize_setup)

Here is a complete example of a soil model devised by Tracy (called a soil context)

Example of the structure used to handle the soil model

/* Parameters defining the unsaturated hydraulic Tracy model */
typedef struct {
double L;
double k_s;
double h_s;
double h_r;
double theta_s;
double theta_r;
/* Array of property values (at cells). These are shared arrays. The
lifecycle of these arrays is not managed by this structure. */
cs_real_t *permeability_values;
cs_real_t *moisture_values;
cs_real_t *capacity_values;
} cs_tracy_param_t;

Example of the first step

cs_real_t bulk_density = 1.0; /* useless here since there is no tracer */
cs_real_t theta_s = 0.45;
cs_gwf_soil_t *s = cs_gwf_add_soil("cells",
bulk_density,
theta_s,
CS_GWF_SOIL_USER);
/* 2.a Create and define a structure to mamange the soil */
cs_tracy_param_t *tp = NULL;
BFT_MALLOC(tp, 1, cs_tracy_param_t);
tp->L = 200;
tp->k_s = 1.15741e-4;
tp->h_s = 0.;
tp->h_r = -100;
tp->theta_s = theta_s;
tp->theta_r = 0.15;
/* Array are set to null since at this stage mesh has not been loaded */
tp->permeability_values = NULL;
tp->moisture_values = NULL;
tp->capacity_values = NULL;
/* 2.b Associate the context structure and the user-defined function
* to the soil */
cs_gwf_soil_set_user(s, /* soil structure */
tp, /* soil context structure */
tracy_update, /* function to update the soil */
tracy_free_context); /* function to free the context */

with the two requested functions (defined for instance as a static function in the file cs_user_parameters.c). These functions have to fullfill the prototype defined in cs_gwf_soil_update_t (for the update of the soil properties) and in cs_gwf_soil_free_context_t (for the free of the soil context).

First, an example of how to update soil properties (tracy_update)

static void
tracy_update(const cs_real_t t_eval,
const cs_mesh_t *mesh,
const cs_cdo_connect_t *connect,
const cs_cdo_quantities_t *quant,
const cs_zone_t *zone,
void *soil_context)
{
CS_UNUSED(mesh);
CS_UNUSED(connect);
CS_UNUSED(quant);
CS_UNUSED(t_eval);
const cs_tracy_param_t *sc = (cs_tracy_param_t *)soil_context;
cs_real_t *head_values = cs_gwf_get_uspf_head_in_law();
const double delta_moisture = sc->theta_s - sc->theta_r;
for (cs_lnum_t i = 0; i < zone->n_elts; i++) {
const cs_lnum_t c_id = zone->elt_ids[i];
const cs_real_t h = head_values[c_id];
const cs_real_t k_r = (h - sc->h_r)/(sc->h_s - sc->h_r);
/* Set the permeability value to the saturated values */
sc->permeability_values[c_id] = sc->k_s * k_r;
/* Set the moisture content (Se = 1 in this case)*/
sc->moisture_values[c_id] = sc->theta_r + k_r * delta_moisture;
/* Set the capacity values */
sc->capacity_values[c_id] = delta_moisture /(sc->h_s - sc->h_r);
} /* Loop on selected cells */
}

and an example of how to free the soil context (tracy_free_context)

static void
tracy_free_context(void **p_soil_context)
{
cs_tracy_param_t *sc = (cs_tracy_param_t *)(*p_soil_context);
BFT_FREE(sc);
*p_soil_context = NULL;
}

Example of the second step

/* 1. Retrieve the soil by its name and then its context */
cs_gwf_soil_t *soil = cs_gwf_soil_by_name("cells");
cs_tracy_param_t *tp = soil->context;
/* 2. Retrieve field values associated to properties to update */
cs_field_t *f = cs_field_by_name("permeability");
assert(f != NULL);
tp->permeability_values = f->val;
/* liquid saturation is also called moisture content */
f = cs_field_by_name("liquid_saturation");
assert(f != NULL);
tp->moisture_values = f->val;
f = cs_field_by_name_try("soil_capacity");
assert(f != NULL);
tp->capacity_values = f->val;

Tracers

The third step (which is not mandatory) is to add tracer(s) thanks to the function cs_gwf_add_tracer This tracer will be advected by the Darcy flux arising from the Richards equation.

There are currently two models :

The first parameter in cs_gwf_add_tracer is a flag which can be built with

Predefined tracers

Here is a simple example for a standard tracer which can be added in the function cs_user_model

{
/*
Add a tracer equation which is unsteady and convected by the darcean flux
This implies the creation of a new equation called eqname along with a
new field called varname.
*/
cs_gwf_tracer_model_t model = 0; /* default model without precipitation
effect */
cs_gwf_tracer_t *tr = cs_gwf_add_tracer(model, /* tracer model */
"Tracer", /* eq. name */
"C"); /* var. name */
/* For "simple" tracer definitions, definition can be made here.
*
* The parameters defining the tracer behavior can be set soil by soil
* (give the soil name as the second argument) or to all soils in one call
* (give NULL as the second argument)
*/
cs_gwf_set_main_tracer_param(tr,
NULL, /* soil name or NULL for all */
0., /* water molecular diffusivity */
1., 0., /* alpha (longi. and transvesal) */
1e-4, /* distribution coef. */
0.01); /* 1st order decay coef. */
}

Remark: Get a tracer structure.

{
/* If the tracer structure has been added but not already defined, one can
* retrieve it thanks to \ref cs_gwf_tracer_by_name
*
* The name of the tracer is the same as the name of the associated
* equation given at the creation of the tracer
*/
cs_gwf_tracer_t *tr = cs_gwf_tracer_by_name("Tracer_01");
}

User-defined tracers

TODO

Automatic postprocessings

It is possible to activate an automatic post-processing of several postprocessings of fields or properties. Here are available flags to activate through the usage of cs_gwf_set_post_options

/* Specify post-processing options */
cs_gwf_set_post_options(CS_GWF_POST_PERMEABILITY |
CS_GWF_POST_DARCY_FLUX_BALANCE |
CS_GWF_POST_DARCY_FLUX_DIVERGENCE |
CS_GWF_POST_DARCY_FLUX_AT_BOUNDARY,
false); /* No reset */

Helper functions

Helper functions for soils

Get a soil structure from its name.

{
/* If a soil structure has been added but not already defined, one can
* retrieve it thanks to \ref cs_gwf_soil_by_name
*
* The name of the soil is the same as the name of the volume zone used at
* the creation of the soil
*/
cs_gwf_soil_t *s1 = cs_gwf_soil_by_name("soil1");
}

There is a similar which retrieve the soil structure from its id (see cs_gwf_soil_by_id).

Helper functions for tracers

Get a tracer structure from its name.

{
/* If the tracer structure has been added but not already defined, one can
* retrieve it thanks to \ref cs_gwf_tracer_by_name
*
* The name of the tracer is the same as the name of the associated
* equation given at the creation of the tracer
*/
cs_gwf_tracer_t *tr = cs_gwf_tracer_by_name("Tracer_01");
}