[FLASH-USERS] reaction-diffusion problem with FLASH?

[Klaus Weide]

On Mon, 13 Nov 2017, Peter Woitke wrote:

> Dear experts,
> 
> is it possible to apply FLASH to a 1D reaction-diffusion problem as
> 
>    d/dt (rho x) = d/dz (rho D dx/dz) + R(x,z)     [1/cm3/s] ?
> 
> The equation describes chemicals in a planetary atmosphere.
> z is the height in the atmosphere [cm], rho the mass density [g/cm3],
> D the diffusion coefficient [m^2/s], and the vector x are the mass
> concentrations of certain gas molecules and solid particles
> [rho x has units 1/cm3]. There is no bulk motion: v=0.

Peter,

You can use FLASH, with some modifications of your own, to set up such a 
problem.

See in particular the documentation of subroutine Grid_advanceDiffusion.
Excerpt:

 DESCRIPTION
      This routine advances a generalized diffusion operator of the form

         A*(df/dt) + C*f = div(B*grad(f)) + D ,

      where
         f = f(x,t) is the  Variable to be diffused (x=1D..3D position);

         A,B,C,D are optional given scalar factors/terms that may depend
         on position; they are either physically constant in time, or at
         least considered time-independent for the purpose of the operation
         implemented here (typically by computing their values from the
         solution state reached by the previous time step).


If you want to try this, I would recommend trying the implementation in 
source/Grid/GridSolvers/HYPRE (using the HYPRE library, obviously) first.

> Both the diffusion constant D(z) and the mass density rho(z) vary
> by about 8 orders of magnitude from the bottom to the top of the atmosphere.
> That's why I need the more general formulation of
> diffusion, can't use  rho D d2x/dz2.

As quoted above, your factor D (called B above) can be location-dependent
for the Grid_advanceDiffusion interface.

You'll have to investigate whether the solver can deal well with your 
requirements. Note that with HYPRE you have a variety of choices for 
solver and preconditioner methods (for the underlying time-independent 
Helmholtz problem that is fed to the HYPRE library to solve).

> There will be a certain fixed temperature structure T(z) assumed,
> and D(z) will be adjusted to turbulent mixing.
> 
> The reaction rates R depend on concentrations x. They are often
> huge compared to the influence of diffusion, hence it's a stiff
> reaction-diffusion problem, and the local chemistry is close to kinetic
> equilibrium. I am interested in how diffusion changes the kinetic equilibrium
> over long timescales.
> 
> I am actually most interested in the t -> infinity solution
> where one could drop the d/dt(rho x) term.

The situation where the solution changes only slowly wrt physical t should
suit the Grid_advanceDiffusion interface.  (Our implementations are not 
really implicit wrt our D term [part of _your_ R term], anyway.)

> I noticed that there is the flash4.4/source/physics/Diffuse
> directory. Does anyone has experience with using that?
> Is there a test problem with diffusion that I could study
> and visualise?

Under source/physics/Diffuse/DiffuseMain you will find various wrappers 
around Grid_advanceDiffusion (ignore for now subdirectories other than 
"Unsplit"). They prepare variables etc. to ultimately call 
Grid_advanceDiffusion to solve a physical diffusive (sub)problem.
(A much more elaborate wrapper is in 
source/physics/RadTrans/RadTransMain/MGD/RadTrans.F90).
None of the existing wrappers there are aware of chemical species, 
composition, or abundances; you may want to construct one, perhaps
starting with a much simplified version of
source/physics/Diffuse/DiffuseMain/Unsplit/diff_advanceTherm.F90
renamed to diff_advanceSpecies ...

Perhaps the simplest test problem is "ConductionDeltaSaDiff" (a 
"standalone-diffusion" [i.e., not flux-based] version of ConductionDelta), 
briefly described in the Users Guide.


> Any advice would be highly appreciated.

HTH,

Klaus