Data setting for the groundwater flow module

Introduction

The Hydrogeology module of \CS is a numerical model for water flow and solute transport in continuous porous media. The flow part is based on the Richards equation, derived from the Darcy law and the conservation of mass. The transport part is based on the the classical advection-diffusion equation, slightly modified to account the specificities of underground transport.

This module can be used to simulate transfers of water and solutes in several saturated and/or unsaturated porous media. The flow part can be steady or unsteady, with scalar or tensorial permeabilities and allows any type of soil water retention model, such as the van Genuchten model. The transport part considers dispersion, sorption and radioactive decay. The partition between soil and water phases can be modeled by a classical Kd approach or an alternative EK (equilibrium-kinetic) model. Additionaly solutes precipitation/dissolution phenomena can also be taken into account by an instantaneous model.

Physical concepts and equations are presented in the theory guide.

The groundwater flow module is recent, and thus, has few limitations:

• The weighted gradient computation, required for tetrahedral meshes and high permeability ratio, can only be used for isotropic soils.
• Only one solute with anisotropic dispersion can be treated.

Activation of the module

The module can be activated in the usppmo routine in cs_user_parameters.f90. The corresponding keyword is idarcy in the Module for calculation options module:

! Set groundwater flow module active (1: active, 0: inactive)
ippmod(idarcy) = 1

Note that the activation of the module requires to desactivation the turbulent model in usipph routine in cs_user_parameters.f90 file:

if (ixmlpu.eq.0) then
  ! For groundwater flow module, turbulent model is set to 0 for security reason
  iturb = 0
endif

Specific parameters

When the module is activated, its specific input parameters should be set in the user_darcy_ini1 routine of cs_user_parameters.f90 file. An example is given in cs_user_parameters-richards.f90.

Flow part

The permeability can be isotropic (scalar) or anisotropic (tensor) but all soils will be treated in the same way (isotropic or anisotropic):

! Set permeability to isotropic (0) or anisotropic (1) for all soils
darcy_anisotropic_permeability = 0

The primary variable of the groundwater flow module is the hydraulic head H=h+z. In order to switch easily to the pressure head, the keyword darcy_gravity can be used and the value darcy_gravity_x/y/z will defined the direction:

! Set gravity to pass from H to h. Example for H = h + z:
darcy_gravity = 1

The convergence criteron of the Newton scheme can be set over pressure or over velocity. It is recommended to keep the criteron over pressure:

! Set convergence criteron of the Newton scheme over pressure (0) or over velocity (1).
! It is recommended to keep the criteron over pressure.
darcy_convergence_criterion = 0

Transport part

The dispersion can be isotropic (scalar) or anisotropic (tensor) but all solutes will be treated in the same way (isotropic or anisotropic):

! Set dispersion to isotropic (0) or anisotropic (1) for all solutes
darcy_anisotropic_dispersion = 0

The transient transport can be based on a steady or unsteasy darcian velocity field:

! Set if the transport part is based on a steady (0) or unsteady (1) flow field
darcy_unsteady = 0

The partition between solid and liquid phases can be modelled by a classical Kd approach or an alternative EK (equilibrium-kinetic) model. Additionally solutes precipitation/dissolution phenomena can also be taken into account by an instantaneous model.

! Get soil-water partition structure.
call field_get_key_struct_gwf_soilwater_partition(ivarfl(isca(1)), &
                                                  sorption_sca1)
 
! Set the sorption model to Kd approach (0) or EK model (1),
! Kd approach is set by default.
sorption_sca1%kinetic = 1
 
! Enable precipitation model, by default, there is no precipitation.
sorption_sca1%imxsol = 0 ! imxsol will hold the solubility index field id
 
! Set the modifications in the soil-water partition structure.
call field_set_key_struct_gwf_soilwater_partition(ivarfl(isca(1)), &
                                                  sorption_sca1)

The radioactive decay is treated as a source term in the transport equation (in covofi.f90). It is set in cs_user_parameters.f90 as follows:

! In groundwater module flow, the radioactive decay of solute is treated as a
! source term in the transport equation
 
! Get radioactive decay rate key
call field_get_key_id("fo_decay_rate", key_decay)
do ii = 1, nscal
  decay_rate = 3.d-1
  ! Set radioactive decay rate
  call field_set_key_double(ivarfl(isca(ii)), key_decay, decay_rate)
enddo

Numerical parameters

Specific numerical parameters can be set in usipsu routine of cs_user_parameters.f90 file. An example is given in cs_user_parameters-richards.f90:

Flow part

In the case of soils of very different permeabilities, the resolution of the Richards equation requires a weighted gradient computation for tetrahedral meshes. This option is only available for soils with isotropic permeabilities (darcy_anisotropic_permeability = 0) for now.

 
call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt)
 
! Set gradient computation to weighted (1) for high permeability ratio in
! tetrahedral meshes.
! Only works with isotropic permeability and the
! standard least square gradient reconstruction.
if (darcy_anisotropic_permeability.eq.0) then
  vcopt%iwgrec = 1
  imrgra = 1
  call field_set_key_struct_var_cal_opt(ivarfl(ipr), vcopt)
endif

It is recommended to choose low criteria for gradient reconstruction in order to obtain a smooth darcian velocity field for the transport part. For instance:

! Set low criteria for gradient reconstruction to obtain a
! smooth velocity field
 
call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt)
 
vcopt%nswrsm = 5000
vcopt%epsrsm = 1.d-10
vcopt%nswrgr = 5000
vcopt%epsrgr = 1.d-10
vcopt%epsilo = 1.d-13
 
call field_set_key_struct_var_cal_opt(ivarfl(ipr), vcopt)
 

Transport part

In the case of soils of very different diffusion (dispersion or molecular diffusion), the resolution of the transport equation requires a weighted gradient computation for tetrahedral meshes.

if (darcy_anisotropic_dispersion.eq.0) then
  do ii = 1, nscal
    ! Set gradient computation to weighted (1) for high permeability ratio in tetrahedral meshes.
    ! Only works with isotropic diffusion.
    call field_get_key_struct_var_cal_opt(ivarfl(isca(ii)), vcopt)
    vcopt%iwgrec = 1
    ! Set higher maximum number of iteration for reconstruction to increase the accuracy
    vcopt%nswrsm = 10
    call field_set_key_struct_var_cal_opt(ivarfl(isca(ii)), vcopt)
  enddo
endif

Time parameters

The total number of iterations and the reference time step are also set in this routine.

ntmabs = 500
dtref  = 1.d-3

However, the time step can be modified in cs_user_extra_operations.f90 (see cs_user_extra_operations-flux.f90) in order to modif the time step with time:

! Example of time modification
do iel = 1, ncel
  dt(iel) = (ttcabs**0.95d0) * 5.d-2
  if (dt(iel).lt.5.d-2) dt(iel) = 5.d-2
  if (dt(iel).gt.1.d3) dt(iel) = 1.d3
enddo

Physical parameters

Physical parameters can be set in usphyv routine of cs_user_physical_properties.f90 file. This section presents two examples that can be found in in cs_user_physical_properties-richards_sat.f90 and cs_user_physical_properties-richards_unsat.f90.

Note that, in the flow part, depending on the variable darcy_anisotropic_permeability, the permeability storage table is permeability (isotropic case) or tensor_permeability (anisotropic case). For the transport part, the isotropic part of the diffusion (molecular diffusion and isotropic dispersion) is always stored in cpro_vscalt. Only anisotropic dispersion (i.e. darcy_anisotropic_diffusion = 1) is stored in the tensor visten.

Case of two saturated soils

This example shows how to set physical parameters for two fully saturated soils (isotropic or anisotropic permeability) and several solutes with isotropic dispersion:

Flow part

! Set parameters for soil 1
call getcel('SOIL1', ncelt, lstcel)
 
! Loop on cell of the list
do icelt = 1, ncelt
 
  ! Cell number
  iel = lstcel(icelt)
 
  ! Set saturation and capacity (0 if saturated)
  saturation(iel) = 0.6d0
  capacity(iel) = 0.d0
 
  ! Set permeability
  if (darcy_anisotropic_permeability.eq.0) then
    permeability(iel) = 1.d-1
  else
    tensor_permeability(1,iel) = 1.d-1
    tensor_permeability(2,iel) = 1.d-1
    tensor_permeability(3,iel) = 1.d-2
    tensor_permeability(4:6,iel) = 0.d0
  endif
 
enddo
 
! Set parameters for soil 2
call getcel('SOIL2', ncelt, lstcel)
 
! Loop on cell of the list
do icelt = 1, ncelt
 
  ! Cell number
  iel = lstcel(icelt)
 
  ! Set saturation and capacity (0 if saturated)
  saturation(iel) = 0.4d0
  capacity(iel) = 0.d0
 
  ! Set permeability
  if (darcy_anisotropic_permeability.eq.0) then
    permeability(iel) = 5.d-1
  else
    tensor_permeability(1,iel) = 5.d-1
    tensor_permeability(2,iel) = 5.d-1
    tensor_permeability(3,iel) = 5.d-2
    tensor_permeability(4:6,iel) = 0.d0
  endif
 
enddo

Transport part

For every solute, the isotropic dispersion should be set in all soils. For instance:

  ! Definition of the isotropic diffusion (dispersion and moleculer diffusion)
  call getcel ('SOIL1', ncelt, lstcel)
  darcy_isotropic_dispersion = 1.d0
  molecular_diffusion = 1.d-6
  do icelt = 1, ncelt
    iel = lstcel(icelt)
    if (ifcvsl.ge.0) then
      velocity_norm = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
      cpro_vscalt(iel) = darcy_isotropic_dispersion * velocity_norm + &
                         saturation(iel) * molecular_diffusion
    endif
  enddo
 
  ! Definition of the isotropic diffusion (dispersion and moleculer diffusion)
  call getcel ('SOIL2', ncelt, lstcel)
  darcy_isotropic_dispersion = 0.2d0
  molecular_diffusion = 1.d-8
  do icelt = 1, ncelt
    iel = lstcel(icelt)
    if (ifcvsl.ge.0) then
      velocity_norm = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
      cpro_vscalt(iel) = darcy_isotropic_dispersion * velocity_norm + &
                         saturation(iel) * molecular_diffusion
    endif
  enddo
 

Unsaturated media

This example shows how to set physical parameters for a single variably saturated soil (isotropic or anisotropic permeability) and a single solute with molecular diffusion, anisotropic dispersivity and sorption. The van Genuchten model, coupled with the Mualem condition, is used to determine the relation between the moisture content and the presure head (h).

Flow part

First the permeability and the van Genuchten parameters are set:

! Set intrinsic permeability (only depends on soil)
if (darcy_anisotropic_permeability.eq.0) then
  ki = 1.d0
else
  ki_xx = 1.d0
  ki_yy = 1.d0
  ki_zz = 1.d-1
endif
 
! Set values of the Van Genuchten model parameters
ks_param = 0.3d0
thetar_param = 0.078d0
thetas_param = 0.3d0
n_param = 1.56d0
m_param = 1-1/n_param !(Mualem condition)
l_param = 0.5d0
alpha_param = 0.036d0

As the van Genuchten law is based on the pressure head (h), the gravity term is added if necessary:

  ! Switch from hydraulic head (H=h+z) to pressure head (h)
  darcy_h = cvar_pr(iel)
  if (darcy_gravity.eq.1) then
    darcy_h = cvar_pr(iel) - xyzcen(3,iel)
  endif

In order to compute the capacity, the saturation and the permeability, the saturated and the unsaturated parts are treated differently.

In the saturated part (h>=0), the water content is equal to the saturated water content, the permeability is equal to the saturated permeability and the capacity is equal to zero:

  ! Saturated part (h<=0)
  if (darcy_h.ge.0) then
 
    capacity(iel) = 0.d0
    saturation(iel) = thetas_param
 
    if (darcy_anisotropic_permeability.eq.0) then
      permeability(iel) = ks_param * ki
    else
      tensor_permeability(1,iel) = ks_param*ki_xx
      tensor_permeability(2,iel) = ks_param*ki_yy
      tensor_permeability(3,iel) = ks_param*ki_zz
      tensor_permeability(4:6,iel) = 0.d0
    endif

In the unsaturated part (h<0), the water content, the capacity and the permeability are function of the pressure head. They are determined thanks to the van Genuchten law:

  ! Unsaturated part (h<0)
  else
 
    tmp_1 = abs(alpha_param*darcy_h)**n_param
    tmp_2 = 1.d0 / (1.d0 + tmp_1)
    se_param = tmp_2**m_param
 
    capacity(iel) = -m_param*n_param*(tmp_1)* &
      (thetas_param-thetar_param)*se_param*tmp_2/darcy_h
    saturation(iel) = thetar_param + se_param*(thetas_param-thetar_param)
 
    if (darcy_anisotropic_permeability.eq.0) then
      permeability(iel) = ks_param*ki*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
    else
      tensor_permeability(1,iel) = ks_param*ki_xx*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
      tensor_permeability(2,iel) = ks_param*ki_yy*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
      tensor_permeability(3,iel) = ks_param*ki_zz*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
      tensor_permeability(4:6,iel) = 0.d0
    endif
 
  endif

Transport part

First, the values of the longitudinal and transversal dispersivity as well as the molecular diffusion of the single solute are set as follow:

! Set values of the longitudinal and transversal dirpersivity
darcy_anisotropic_dispersion_l = 2.d0
darcy_anisotropic_dispersion_t = 1.d-1
 
! Set value of the molecular diffusion
molecular_diffusion = 1.d-3

The molecular diffusion (isotropic term) is stored in cpro_vscalt and computed as:

! Computation of molecular diffusion of the diffusion term
if (ifcvsl.ge.0) then
  call field_get_val_s(ifcvsl, cpro_vscalt)
  do iel = 1, ncel
    cpro_vscalt(iel) = saturation(iel)*molecular_diffusion
  enddo
else
  cpro_vscalt => null()
endif

The anisotropic dispersion is stored in visten and computed as:

! Computation of the isotropic dispersivity
do iel = 1, ncel
 
  ! Computation of the norm of the velocity
  velocity_norm = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
 
  ! Tensorial dispersion is stored in visten
  tmp_lt = darcy_anisotropic_dispersion_l-darcy_anisotropic_dispersion_t
  visten(1,iel) = darcy_anisotropic_dispersion_t*velocity_norm + tmp_lt*vel(1,iel)**2/(velocity_norm+1.d-15)
  visten(2,iel) = darcy_anisotropic_dispersion_t*velocity_norm + tmp_lt*vel(2,iel)**2/(velocity_norm+1.d-15)
  visten(3,iel) = darcy_anisotropic_dispersion_t*velocity_norm + tmp_lt*vel(3,iel)**2/(velocity_norm+1.d-15)
  visten(4,iel) = tmp_lt*vel(2,iel)*vel(1,iel)/(velocity_norm+1.d-15)
  visten(5,iel) = tmp_lt*vel(2,iel)*vel(3,iel)/(velocity_norm+1.d-15)
  visten(6,iel) = tmp_lt*vel(3,iel)*vel(1,iel)/(velocity_norm+1.d-15)
 
enddo

Sorption parameters like Kd, kplus or kminus (in case of EK model) can be set as follows. Soil density is also needed to compute the transport delay. If precipitation is taken into account, the solubility index needs to be set as well.

! Set soil density (bulk density!) for delay computation (delay = 1 + soil_density * K_d / saturation)
do iel = 1, ncel
  soil_density(iel) = 1.5d0
enddo
 
! Get soil-water partition structure
call field_get_key_struct_gwf_soilwater_partition(ivarfl(isca(1)), &
                                                  sorption_scal)
 
! Index field for kd
call field_get_val_s(sorption_scal%ikd, cpro_kd)
 
! Index field for EK model parameters (kplus and kminus)
call field_get_val_s(sorption_scal%ikp, cpro_kplus)
call field_get_val_s(sorption_scal%ikm, cpro_kminus)
 
! Set sorption parameters
do iel=1, ncel
  cpro_kd(iel) = 5.d0
  ! if EK model is chosen, set specific parameters
  cpro_kplus(iel) =  1.d-3
  cpro_kminus(iel) = 1.d-4
enddo
 
!Index field for cpro_mxsol index (if precipitation option is activated)
call field_get_val_s(sorption_scal%imxsol, cpro_mxsol)
 
do iel=1, ncel
  cpro_mxsol(iel) = 10.d0
enddo
 

Source terms

Source terms can be added to scalar transport equation in cs_user_source_terms function of the cs_user_source_terms.c file.

Chemicals release

Substances can be gradually released within the soil.

To define progressive leaching, the following local variables initialization can be used:

  /* local variables */
 
  const cs_field_t  *f = cs_field_by_id(f_id);
  const cs_real_t  t = cs_glob_time_step->t_cur;
 
  const cs_real_t  *cell_f_vol = cs_glob_mesh_quantities->cell_vol;
 
  /* Map saturation and delay */
 
  char f_name[64];
  snprintf(f_name, 63, "%s_delay", f->name); f_name[63] = '\0';
  const cs_real_t  *delay = cs_field_by_name(f_name)->val;
 
  const cs_real_t  *saturation = cs_field_by_name("saturation")->val;

Additional logging based on the field's general verbosity can be done in the usual manner if desired:

  const cs_var_cal_opt_t  *var_cal_opt
    = cs_field_get_key_struct_const_ptr(f,
                                        cs_field_key_id("var_cal_opt"));
  if (var_cal_opt->iwarni >= 1)
    bft_printf(" User source terms for variable %s\n",
               cs_field_get_label(f));

Then, define the source terms for leaching from a volume zone named "LEACHING_ZONE":

  cs_real_t lambda = 1.e-2;        /* first order decay coefficient */
  cs_real_t leaching_time = 1.e2;  /* leaching duration */
 
  const cs_zone_t *z = cs_volume_zone_by_name_try("LEACHING_ZONE");
 
  /* progressive leaching */
 
  if (z != NULL && t < leaching_time) {
    cs_real_t leaching_volume = z->f_measure;
 
    for (cs_lnum_t i = 0; i < z->n_elts; i++) {
      cs_lnum_t c_id = z->elt_ids[i];
 
      st_exp[c_id] =   cell_f_vol[i] * exp(-lambda*t)
                     / (  leaching_time * leaching_volume
                        * saturation[i] * delay[i]);
    }
  }

Initialization

The initialization of the variables required for the flow part (hydraulic head H) and transport part (concentration c) can be done globally:

! Global initialisation
  do iel = 1, ncel
 
!   Initialisation of the hydraulic pressure H with a contant gradient of 1
!   among z axis and -0.01 among x axis
    cvar_pr(iel) = 1.d0*xyzcen(3,iel) - 1.d-2*xyzcen(1,iel)
 
!   Null concentration by default
    cvar_scal_1(iel) = 0.d0
 
!   Velocity initialisation for security reason
    cvar_vel(1,iel) = 0.d0
    cvar_vel(2,iel) = 0.d0
    cvar_vel(3,iel) = 0.d0
  enddo

or by selecting a precise soil:

! Initialisation per group of volume
  call getcel('SOURCE', ncelt, lstelt)
  do icelt = 1, ncelt
    iel = lstelt(icelt)
    cvar_scal_1(iel) = 1.d0
  enddo

If EK model is considered, sorbed concentration must be initialized.

  ! Index field of sorbed concentration
  call field_get_key_id("gwf_sorbed_concentration_id", keysrb)
  call field_get_key_int(ivarfl(isca(1)), keysrb, isorb)
  call field_get_val_s(isorb, cpro_sorb)
 
  do iel = 1, ncel
    ! no initial contamination of sorbed phase
    cpro_sorb(iel) = 0.d0
  enddo

If precipitation phenomenon is taken into account, precipitated concentration has to be initialized.

  ! Index field of precipitated concentration
  call field_get_key_id("gwf_precip_concentration_id", keypre)
  call field_get_key_int(ivarfl(isca(1)), keypre, igwfpr)
  call field_get_val_s(igwfpr, cpro_precip)
 
  do iel = 1, ncel
    ! no initial precipitation phase
    cpro_precip(iel) = 0.d0
  enddo

Boundary conditions

For groundwater flows of water and solutes, the undefined type face iindef is used to impose Dirichlet, Neumann and mixed boundary conditions on hydraulic head H (here pressure) and solutes. Several examples can be found in cs_user_boundary_conditions-richards.f90.

Dirichlet boundary conditions

Dirichlet boundary conditions can be used to impose a value for the hydraulic head H and the concentration c at a given boundary face:

call getfbr('FACE1', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
! Face number
  ifac = lstelt(ilelt)
 
! Undefined type face
  itypfb(ifac) = iindef
 
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
  icodcl(ifac,iu) = 3
  icodcl(ifac,iv) = 3
  icodcl(ifac,iw) = 3
  rcodcl(ifac,iu,3) = 0.d0
  rcodcl(ifac,iv,3) = 0.d0
  rcodcl(ifac,iw,3) = 0.d0
 
! Dirichlet BC on hydraulic head (H = h + z) to impose a constant value
  icodcl(ifac,ipr) = 1
  rcodcl(ifac,ipr,1) = 10.d0
 
! Dirichlet BC on centration C (g/m^3)
  do ii = 1, nscal
    icodcl(ifac,isca(ii)) = 1
    rcodcl(ifac,isca(ii),1) = 1.d0
  enddo
enddo

It can also be used to impose a hydraulic head profile at another face:

call getfbr('FACE2', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
! Face number
  ifac = lstelt(ilelt)
 
! Undefined type face
  itypfb(ifac) = iindef
 
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
  icodcl(ifac,iu) = 3
  icodcl(ifac,iv) = 3
  icodcl(ifac,iw) = 3
  rcodcl(ifac,iu,3) = 0.d0
  rcodcl(ifac,iv,3) = 0.d0
  rcodcl(ifac,iw,3) = 0.d0
 
! Dirichlet BC on hydraulic head (H = h + z) to impose a constant gradient over z axis
! here \f$ \grad H \cdot \vect{z} = -0.1 \f$
  icodcl(ifac,ipr) = 1
  rcodcl(ifac,ipr,1) = -1.d0 * cdgfbo(3,ifac)
enddo

Neumann boundary conditions

Neumann boundary conditions can be used to impose fluxes at boundaries:

call getfbr('FACE3', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
! Face number
  ifac = lstelt(ilelt)
 
! Undefined type face
  itypfb(ifac) = iindef
 
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
  icodcl(ifac,iu) = 3
  icodcl(ifac,iv) = 3
  icodcl(ifac,iw) = 3
  rcodcl(ifac,iu,3) = 0.d0
  rcodcl(ifac,iv,3) = 0.d0
  rcodcl(ifac,iw,3) = 0.d0
 
! Neumann BC on hydraulic head (H = h + z) to impose a gradient among the surface normal
! Here \f$ \grad H \cdot \vect{n} = 0 \f$
  icodcl(ifac,ipr) = 3
  rcodcl(ifac,ipr,3) = 0.d0
 
! Neumann BC on centration C for boundary surface with outward or null normal flow $V_out$:
! It allows to impose a gradient among the surface normal for a diffusive flux as the
! convective flux is defined by the $C*V_out$.
! Here \f$ \grad C \cdot \vect{n} = 0 \f$
  do ii = 1, nscal
    icodcl(ifac,isca(ii)) = 3
    rcodcl(ifac,isca(ii),3) = 0.d0
  enddo
enddo

Note that, for the transport part, Neumann boundary conditions can only be used for boundary surface with outward or null normal flow. In both cases, the prescribed flux is the diffusive flux.

Mixed boundary conditions

The mixed boundary conditions (Robin) can be used to impose a concentration flux at an entrance (inward normal flow at a boundary face). The following example explains how to determine the two parameters of the mixed boundary in order to impose a total flux:

call getfbr('FACE4', nlelt, lstelt)
 
! Pointers to the mass flux
call field_get_key_int(ivarfl(ipr), kbmasf, iflmab)
call field_get_val_s(iflmab, bmasfl)
 
do ilelt = 1, nlelt
 
! Face number
  ifac = lstelt(ilelt)
 
! Undefined type face
  itypfb(ifac) = iindef
 
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
  icodcl(ifac,iu) = 3
  icodcl(ifac,iv) = 3
  icodcl(ifac,iw) = 3
  rcodcl(ifac,iu,3) = 0.d0
  rcodcl(ifac,iv,3) = 0.d0
  rcodcl(ifac,iw,3) = 0.d0
 
! Dirichlet BC on hydraulic head (H = h + z) to impose a constant value
  icodcl(ifac,ipr) = 1
  rcodcl(ifac,ipr,1) = 10.d0
 
! In order to impose a radioactive activity R_act (Bq) at a surface S_act (m^2)
! during a period of time T_act, a mixed (or Robin) BC is used.
! To do so, two quantities are required:
! * the velocity at an entrance is Vel = - mass_flux / surf
!   or Vel = mass_flux / surf at an exit
! * the reference concentration Cref = R_act / (S_act * dt * Vel)
  r_act = 0.5d0
  s_act = 50.d0
  t_act = 10.d0
 
! The total flux is imposed from 0 to T_act
  vel = - bmasfl(ifac)/surfbn(ifac)
  cref = r_act / (s_act * dt(ntcabs) * vel)
  if (ttcabs.gt.t_act) cref = 0.d0
  do ii = 1, nscal
    icodcl(ifac,isca(ii)) = 1
    rcodcl(ifac,isca(ii),1) = cref
    rcodcl(ifac,isca(ii),2) = vel
  enddo
enddo

Flux computations

In order to compute fluxes at innner or boundary surfaces, the file cs_user_extra_operations.f90 can be used. An example can be found in cs_user_extra_operations-richards_flux.f90. It shows how to compute the total fluxes at to different surface and how to write the evolution with time:

pas_iter = 10 ! Number of iterations separating two flux calculations.
 
if (mod(ntcabs,pas_iter).eq.0) then
 
  iprev = 0
  inc = 1
  iccocg = 1
  eps_geom = 1.d-10
 
  call field_get_key_int(ivarfl(iu), kimasf, iflmas)
  call field_get_key_int(ivarfl(iu), kbmasf, iflmab)
  call field_get_val_s(iflmas, imasfl)
  call field_get_val_s(iflmab, bmasfl)
  call field_get_key_int(ivarfl(isca(1)), kivisl, ifcvsl)
  call field_get_val_s(ifcvsl, cpro_vscalt)
  scalar_diffusion(1:ncel) = cpro_vscalt(1:ncel)
  call viscfa (imvisf, scalar_diffusion, viscf, viscb)
 
! Example of tracer flux computation through an internal surface.
! Fluxes are calculated whithout reconstruction, and with a simple definition of the concentration at faces
! (i.e. mean of the two adjacent concentrations).
! Thus, the fluxes are not coherent with the real fluxes calculated in bilsc2 (variable "flux")
! They are approximated and should not be used to verify the exact conservation of mass
 
  call getfac('PLANE1', nlelt, lstelt)
 
  flux_in = 0.d0
  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    iel11 = ifacel(1,ifac)
    iel22 = ifacel(2,ifac)
 
!   Test locality of the cell (only for mpi computation)
    if (iel22.le.ncel) then
      tra_face = 5.d-1*(cvar_scal_1(iel11)+cvar_scal_1(iel22))
      tra_diff = cvar_scal_1(iel11) - cvar_scal_1(iel22)
    else
!     Remove duplicated boundary boundary faces for parallel computation
      tra_face = 0.d0
      tra_diff = 0.d0
    endif
 
    flux_tmp = imasfl(ifac)*tra_face + viscf(ifac)*tra_diff
 
!   We impose a norm on the direction of the flux, based on the direction
!   of main coordinates, to be sure that fluxes of all faces are all summed
!   in the same direction
    if ( (surfac(1,ifac).le.(-eps_geom) ) .or.   &
       ( (abs(surfac(1,ifac)).le.eps_geom).and.  &
         (surfac(2,ifac).le.(-eps_geom)) ) .or.  &
       ( (abs(surfac(1,ifac)).le.eps_geom).and.  &
         (abs(surfac(2,ifac)).le.eps_geom).and.  &
         (surfac(3,ifac).le.(-eps_geom)) ) ) then
      flux_tmp = -flux_tmp
    endif
    flux_in = flux_in + flux_tmp
  enddo
 
! Example of boundary flux computation
! Fluxes are calculated whithout reconstruction, and with a simple definition of the concentration at faces
! (i.e. mean of the two adjacent concentrations).
 
  call getfac('TOP_SOIL1', nlelt, lstelt)
 
  flux_out = 0.d0
  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    iel11 = ifacel(1,ifac)
    iel22 = ifacel(2,ifac)
 
!   Test locality of the cell
    if (iel22.le.ncel) then
      tra_face = 5.d-1*(cvar_scal_1(iel11)+cvar_scal_1(iel22))
      tra_diff = cvar_scal_1(iel11) - cvar_scal_1(iel22)
    else
!     Remove duplicated boundary boundary faces for parallel computation
      tra_face = 0.d0
      tra_diff = 0.d0
    endif
!
    flux_tmp = imasfl(ifac)*tra_face + viscf(ifac)*tra_diff
 
!   We impose a norm on the direction of the flux, based on the direction
!   of main coordinates, to be sure that fluxes of all faces are all summed
!   in the same direction
    if ( (surfac(1,ifac).le.(-eps_geom) ) .or.   &
       ( (abs(surfac(1,ifac)).le.eps_geom).and.  &
         (surfac(2,ifac).le.(-eps_geom)) ) .or.  &
       ( (abs(surfac(1,ifac)).le.eps_geom).and.  &
         (abs(surfac(2,ifac)).le.eps_geom).and.  &
         (surfac(3,ifac).le.(-eps_geom)) ) ) then
      flux_tmp = -flux_tmp
    endif
    flux_out = flux_out + flux_tmp
  enddo
 
  ! Compute sum for parallel computations
  if (irangp.ge.0) then
    call parsom(flux_in)
    call parsom(flux_out)
  endif
 
  ! Write fluxes in file
  impout = impusr(2)
  if (irangp.le.0) then
    open(impout,file='flux_tracer1.dat',position="append")
    write(impout,97) ttcabs, flux_in, flux_out
97      format(3g17.9)
    close(impout)
  endif
 
endif