AdvReactDensity.md

Solving the Mass Conservation of Snow/Firn

General Description

This page explains how to use the general AdvectionReactionSolver from the Elmer distribution to get the density evolution in case of a compressible material (snow/firn) under a given velocity field computed from the Porous Solver. The AdvectionReactionSolver solves the general equation

{{\partial A}/{\partial t}} + div (A u) + gamma A=sigma

where u is the velocity vector. In the particular case of the mass conservation equation, one has therefore gamma = 0 and sigma = 0. Solving for the true density (kg/m^3) or the relative density is equivalent (but limit values and Dirichlet boundary conditions have to be set accordingly).

Note 1: equation (4.1) in the Elmer Model Manual for the AdvectionReaction solver is not correct. The previous equation is the one implemented.

Note 2: Have a look to this post on the Elmer Forum regarding the initialisation of both the DG primary and exported variables of the AdvectionReaction solver (see the example at the end of this page).

SIF contents

The Solver section looks like:

Solver 8
  Equation = "AdvReact"
  Exec Solver = "After Timestep"
  Procedure = File "AdvectionReaction" "AdvectionReactionSolver"
  ! this is the DG variable, which is not part of the output
  Variable =  -nooutput "DGdens"
  ! this tells that the solver is run on DG mesh
  Discontinuous Galerkin = Logical True
  ! the solver can account for upper and lower limits of the variable
  ! imposed by formulation of an variational inequality (VI)
  ! next line switches the VI to be accounted for
  Limit Solution = Logical True

  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations  = 1000
  Linear System Preconditioning = ILU1
  Linear System Convergence Tolerance = 1.0e-06
  ! Variational inequality makes it a non-linear problem
  Nonlinear System Max Iterations = 40
  Nonlinear System Min Iterations = 2
  Nonlinear System Convergence Tolerance = 1.0e-04

  ! This is the variable that is used to interpolate
  ! the DG solution to the regular FEM mesh in order
  ! to get a correct output
  Exported Variable 1 = Relative Density
  Exported Variable 1 DOFS = 1
End

The source in case of the mass conservation equation is 0

Body Force 1
 ...
  DGDens Source = Real 0.0 
 End

Initial and boundary conditions are then being set for the DG variable and not the exported one!

Initial Condition 1
  ...
  DGDens = Real 0.4
End

! only Dirichlet BC can be set
! the solver automatically uses this
! condition only on inflow boundaries
! outflow boundaries are ignored
Boundary Condition 2
  Name = "surf"
  Target Boundaries = 2
  Body ID = 2
  ...
  ! relative density on the upper surface
  DGDens = Real 0.4
End

The Material section contains a zero reaction rate as well as the upper/lower limit for the DG variable

Material 1
 ..
 ! Relative density must stay < 1
 DGDens Upper Limit = Real 1.0

 ! a minimum relative density is recommended for the Porous solver 
 DGDens Lower Limit = Real 0.3

 !Reaction rate is equal to zero
 DGDens Gamma = Real 0.0
End

Examples

A 1D example build from an analytical solution can be found in [ELMER_TRUNK]/elmerice/Tests/Density. In that case, the velocity and density are inversely proportional (u(z) = K/D(z)).

A 3D example coupling the Porous Solver and the calculation of the density field can be found in [ELMER_TRUNK]/elmerice/Tests/DGsolver.