Varying mathematical model formulations in space and time: a generalized approach

Abstract

Within a natural system, it is common to find different physical and chemical processes in different parts of the domain. The spatial distribution of these processes may also vary in time. Examples include: groundwater flow problems where a part of the domain becomes unsaturated, and solute transport problems that have nonzero concentrations in only a part of the domain. When numerically modeling such systems with codes not designed for this task, this problem can be dealt with by solving trivial equations in all the domain (eg, solving transport over all the domain, even where it is known that no solute is present) or by solving overly general equations (eg, formulating an unsaturated flow problem even there where it is saturated), usually more expensive to solve. However, a more elegant and efficient approach is to adapt the mathematical model in space and time to the characteristics of the system. We present here a design to accomplish this feat. Specifically, our design is of a group of classes that will operate within Proost, an object oriented framework for numerically solving partial differential equations (PDE0s) for hydrogeological problems, with (for now) the finite element method. Both the solution variables of these PDE’s and the coefficients that appear in them are scalar, vector or tensor fields. In addition to the PDE’s these fields can be related by local or constitutive laws. For example, groundwater flow is represented mathematically by the flow equation (a PDE), while the relative permeability that appears in it is a local equation.In the abstraction underlying this design, we distinguish four types of equations that need solving. First, equations that have no relation to the system state, for example the intrinsic permeability of the porous medium as a function of space. These can be calculated straightforwardly. Second, the PDE themselves. They need to be solved with a numerical method, and often also a linearization method such as Newton or Picard. Third, local equations (equations that involve only quantities measurable at the same location where the equation is solved, eg, density as a function of concentration). They can be in explicit or implicit form. In the latter case, they are solved simultaneously as a set of equations in residual form. Fourth, nonlocal equations (for example, fluid density at elements as an average of fluid density at nodes). They can only be explicit and, hence, can be substituted straightforwardly into the PDEs. The most challenging part is solving PDE0s for different unknowns in different parts of the domain. Each PDE, when discretized, leads to a mass balance equation involving fluxes between different nodes (eg, node i and node j). However, we may be solving different PDE’s in node i and node j (or the same PDE but for a different unknown). Therefore, it is necessary to express fluxes from node i to node j as a function of the unknowns of both node i and of node j. With this functionality implemented generically, PDE0s can be solved as a function of different unknowns in different parts of the domain. The constitutive laws needed to close the system can be solved by Newton Rhapson (if no explicit functions are available), with a different set of local equations in different parts of the domain. This flexibility will allow enhancing convergence and reduce calculation time.