mod_ard_phys.t

Go to the documentation of this file.

!> Module containing the physics routines for advection-reaction-diffusion equations
!>
!> This module can be seen as an extension of the reaction-diffusion (rd) module
!> and includes the same reaction systems and more: the Gray-Scott model, the
!> Schnakenberg model, the Brusselator model, the diffusive logistic equation,
!> an analytical testcase from "Numerical solution of time-dependent advection-
!> diffusion-reaction equations" by Hundsdorfer & Verwer, the Oregonator model,
!> the extended Brusselator model, the diffusive Lorenz system and the advection-
!> diffusion equation. See the documentation of the advection-reaction-diffusion
!> module for more information.
!>
!> An advection term can be aplied to these systems of the form: 
!> nabla( (A1/adv_pow) * u^(adv_pow) )   (for the first unknown)
!> nabla( (A2/adv_pow) * v^(adv_pow) )   (for the second unknown, if applicable)
!> nabla( (A3/adv_pow) * w^(adv_pow) )   (for the third unknown, if applicable)
!>
!> IMEX methods are also supported. The implicit system is solved by a
!> multigrid solver coupled into MPI-AMRVAC.
!>
module mod_ard_phys
  use mod_multigrid_coupling
  use mod_comm_lib, only: mpistop
 
  implicit none
  private
 
  !> Indices of the unknowns
  integer, protected, public :: u_ = 1
  integer, protected, public :: v_ = 2 !< For 2 or more equations
  integer, protected, public :: w_ = 3 !< For 3 or more equations
 
  !> Whether particles module is added
  logical, public, protected :: ard_particles = .false.
 
  !> Parameter with which to multiply the reaction timestep restriction
  double precision, public, protected :: dtreacpar = 0.5d0
 
  !> Name of the system to be solved
  character(len=20), public, protected :: equation_name = "gray-scott"
  integer            :: number_of_species  = 2
  integer            :: equation_type      = 1
  integer, parameter :: eq_gray_scott      = 1 ! Gray-Scott model
  integer, parameter :: eq_schnakenberg    = 2 ! Schnakenberg model
  integer, parameter :: eq_brusselator     = 3 ! Brusselator model
  integer, parameter :: eq_logistic        = 4 ! Logistic equation
  integer, parameter :: eq_analyt_hunds    = 5
  integer, parameter :: eq_belousov_fn     = 6 ! Field-Noyes model, or Oregonator
  integer, parameter :: eq_ext_brusselator = 7 ! Extended Brusselator
  integer, parameter :: eq_lorenz          = 8 ! Lorenz system
  integer, parameter :: eq_no_reac         = 9 ! Advection-diffusion equation
 
  !> Diffusion coefficient for first species (u)
  double precision, public, protected :: d1 = 0.05d0
  !> Diffusion coefficient for second species (v) (if applicable)
  double precision, public, protected :: d2 = 1.0d0
  !> Diffusion coefficient for third species (w) (if applicable)
  double precision, public, protected :: d3 = 1.0d0
 
  !> Power of the unknown in the advection term (1 for linear)
  integer, public, protected :: adv_pow = 1
 
  !> Advection coefficients for first species (u)
  double precision, public, protected :: a1(^nd) = 0.0d0
  !> Advection coefficients for second species (v) (if applicable)
  double precision, public, protected :: a2(^nd) = 0.0d0
  !> Advection coefficients for third species (w) (if applicable)
  double precision, public, protected :: a3(^nd) = 0.0d0
 
  !> Parameters for Schnakenberg model
  double precision, public, protected :: sb_alpha = 0.1305d0
  double precision, public, protected :: sb_beta  = 0.7695d0
  double precision, public, protected :: sb_kappa = 100.0d0
 
  !> Feed rate for Gray-Scott model
  double precision, public, protected :: gs_f = 0.046d0
  !> Kill rate for Gray-Scott model
  double precision, public, protected :: gs_k = 0.063d0
 
  !> Parameters for Brusselator model
  double precision, public, protected :: br_a = 4.5d0
  double precision, public, protected :: br_b = 8.0d0
  double precision, public, protected :: br_c = 1.0d0
  double precision, public, protected :: br_d = 1.0d0
 
  !> Parameter for logistic model (Fisher / KPP equation)
  double precision, public, protected :: lg_lambda = 1.0d0
 
  !> Parameters for the Field-Noyes model of the Belousov-Zhabotinsky reaction
  double precision, public, protected :: bzfn_epsilon = 1.0d0
  double precision, public, protected :: bzfn_delta   = 1.0d0
  double precision, public, protected :: bzfn_lambda  = 1.0d0
  double precision, public, protected :: bzfn_mu      = 1.0d0
 
  !> Parameter for Lorenz system (Rayleigh number)
  double precision, public, protected :: lor_r     = 28.0d0
  !> Parameter for Lorenz system (Prandtl number)
  double precision, public, protected :: lor_sigma = 10.0d0
  !> Parameter for Lorenz system (aspect ratio of the convection rolls)
  double precision, public, protected :: lor_b     = 8.0d0 / 3.0d0
 
  !> Whether to handle the explicitly handled source term in split fashion
  logical :: ard_source_split = .false.
 
  !> Boundary condition information for the multigrid method
  type(mg_bc_t), public :: ard_mg_bc(3, mg_num_neighbors)
 
  ! Public methods
  public :: ard_phys_init
 
contains
 
  !> Read this module's parameters from a file
  subroutine ard_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /ard_list/ d1, d2, d3, adv_pow, a1, a2, a3, sb_alpha, sb_beta, sb_kappa, &
        gs_f, gs_k, br_a, br_b, br_c, br_d, lg_lambda, bzfn_epsilon, bzfn_delta, &
        bzfn_lambda, bzfn_mu, lor_r, lor_sigma, lor_b, equation_name, ard_particles, &
        ard_source_split, dtreacpar
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status='old')
       read(unitpar, ard_list, end=111)
111    close(unitpar)
    end do
 
    !> Set the equation type and number of species
    select case (equation_name)
    case ("gray-scott")
       equation_type = eq_gray_scott
       number_of_species = 2
    case ("schnakenberg")
       equation_type = eq_schnakenberg
       number_of_species = 2
    case ("brusselator")
       equation_type = eq_brusselator
       number_of_species = 2
    case ("logistic")
       equation_type = eq_logistic
       number_of_species = 1
    case ("analyt_hunds")
       equation_type = eq_analyt_hunds
       number_of_species = 1
    case ("belousov_fieldnoyes")
       equation_type = eq_belousov_fn
       number_of_species = 3
    case ("ext_brusselator")
       equation_type = eq_ext_brusselator
       number_of_species = 3
    case ("lorenz")
       equation_type = eq_lorenz
       number_of_species = 3
    case ("no_reac")
       equation_type = eq_no_reac
       number_of_species = 1
    case default
       call mpistop("Unknown equation_name (valid: gray-scott, schnakenberg, ...)")
    end select
 
  end subroutine ard_params_read
 
  !> Write this modules parameters to a snapshot
  subroutine ard_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer, parameter                  :: n_par = 0
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er
    integer                             :: idim
 
    call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
  end subroutine ard_write_info
 
  subroutine ard_phys_init()
    use mod_global_parameters
    use mod_physics
    use mod_multigrid_coupling
    use mod_particles, only: particles_init
 
    call ard_params_read(par_files)
 
    physics_type = "ard"
    phys_energy  = .false.
    use_particles = ard_particles
 
    allocate(start_indices(number_species),stop_indices(number_species))
    ! set the index of the first flux variable for species 1
    start_indices(1)=1
    ! Use the first variable as a density
    u_ = var_set_fluxvar("u", "u")
    if (number_of_species >= 2) then
        v_ = var_set_fluxvar("v", "v")
    end if
    if (number_of_species >= 3) then
        w_ = var_set_fluxvar("w", "w")
    end if
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
 
    ! Check whether custom flux types have been defined
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nw))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    phys_get_cmax          => ard_get_cmax
    phys_get_cbounds       => ard_get_cbounds
    phys_get_flux          => ard_get_flux
    phys_to_conserved      => ard_to_conserved
    phys_to_primitive      => ard_to_primitive
    phys_add_source_geom   => ard_add_source_geom
    phys_add_source        => ard_add_source
    phys_get_dt            => ard_get_dt
    phys_write_info        => ard_write_info
    phys_check_params      => ard_check_params
    phys_implicit_update   => ard_implicit_update
    phys_evaluate_implicit => ard_evaluate_implicit
 
    ! Initialize particles module
    if (ard_particles) then
       call particles_init()
    end if
 
  end subroutine ard_phys_init
 
  subroutine ard_check_params
    use mod_global_parameters
    integer :: n, i, iw, species_list(number_of_species)
 
    if (use_imex_scheme) then
       use_multigrid = .true.
       select case(number_of_species)
          case(1)
             species_list = [u_]
             if (d1 == 0.0d0) then
                write(*, *) "Diffusion coefficient cannot be zero in IMEX scheme"
                call mpistop("Zero diffusion in IMEX scheme")
             end if
          case(2)
             species_list = [u_, v_]
             if ((d1 == 0.0d0) .or. (d2 == 0.0d0)) then
                write(*, *) "Diffusion coefficient cannot be zero in IMEX scheme"
                call mpistop("Zero diffusion in IMEX scheme")
             end if
          case(3)
             species_list = [u_, v_, w_]
             if ((d1 == 0.0d0) .or. (d2 == 0.0d0) .or. (d3 == 0.0d0)) then
                write(*, *) "Diffusion coefficient cannot be zero in IMEX scheme"
                call mpistop("Zero diffusion in IMEX scheme")
             end if
       end select
 
       do i = 1, number_of_species
          iw = species_list(i)
 
          ! Set boundary conditions for the multigrid solver
          do n = 1, 2*ndim
             select case (typeboundary(iw, n))
             case (bc_symm)
                ! d/dx u = 0
                ard_mg_bc(i, n)%bc_type = mg_bc_neumann
                ard_mg_bc(i, n)%bc_value = 0.0_dp
             case (bc_asymm)
                ! u = 0
                ard_mg_bc(i, n)%bc_type = mg_bc_dirichlet
                ard_mg_bc(i, n)%bc_value = 0.0_dp
             case (bc_cont)
                ! d/dx u = 0
                ard_mg_bc(i, n)%bc_type = mg_bc_neumann
                ard_mg_bc(i, n)%bc_value = 0.0_dp
             case (bc_periodic)
                ! Nothing to do here
             case (bc_special)
                if (.not. associated(ard_mg_bc(i, n)%boundary_cond)) then
                   write(*, "(A,I0,A,I0,A)") "typeboundary(", iw, ",", n, &
                        ") is 'special', but the corresponding method " // &
                        "ard_mg_bc(i, n)%boundary_cond is not set"
                   call mpistop("ard_mg_bc(i, n)%boundary_cond not set")
                end if
             case default
                write(*,*) "ard_check_params warning: unknown boundary type"
                ard_mg_bc(i, n)%bc_type = mg_bc_dirichlet
                ard_mg_bc(i, n)%bc_value = 0.0_dp
             end select
          end do
       end do
    end if
 
  end subroutine ard_check_params
 
  subroutine ard_to_conserved(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
 
    ! Do nothing (primitive and conservative are equal for ard module)
  end subroutine ard_to_conserved
 
  subroutine ard_to_primitive(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
 
    ! Do nothing (primitive and conservative are equal for ard module)
  end subroutine ard_to_primitive
 
  subroutine ard_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
    use mod_global_parameters
    integer, intent(in)                       :: ixi^l, ixo^l, idim
    double precision, intent(in)              :: w(ixi^s, nw), x(ixi^s, 1:^nd)
    double precision, intent(inout)           :: cmax(ixi^s)
 
    cmax(ixo^s) = abs(a1(idim) * w(ixo^s,u_)**(adv_pow-1))
    if (number_of_species >= 2) then
       cmax(ixo^s) = max(cmax(ixo^s), abs(a2(idim) * w(ixo^s,v_)**(adv_pow-1)))
    end if
    if (number_of_species >= 3) then
       cmax(ixo^s) = max(cmax(ixo^s), abs(a3(idim) * w(ixo^s,w_)**(adv_pow-1)))
    end if
 
  end subroutine ard_get_cmax
 
  subroutine ard_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed, cmax, cmin)
    use mod_global_parameters
    use mod_variables
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wlc(ixi^s, nw), wrc(ixi^s,nw)
    double precision, intent(in)    :: wlp(ixi^s, nw), wrp(ixi^s,nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
    double precision, intent(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
 
    double precision :: wmean(ixi^s,nw)
 
    ! Since the advection coefficient can depend on unknowns,
    ! some average over the left and right state should be taken
    wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
 
    if (present(cmin)) then
       cmin(ixo^s,1) = min(a1(idim) * wmean(ixo^s,u_)**(adv_pow-1), zero)
       cmax(ixo^s,1) = max(a1(idim) * wmean(ixo^s,u_)**(adv_pow-1), zero)
       if (number_of_species >= 2) then
          cmin(ixo^s,1) = min(cmin(ixo^s,1), a2(idim) * wmean(ixo^s,v_)**(adv_pow-1))
          cmax(ixo^s,1) = max(cmax(ixo^s,1), a2(idim) * wmean(ixo^s,v_)**(adv_pow-1))
       end if
       if (number_of_species >= 3) then
          cmin(ixo^s,1) = min(cmin(ixo^s,1), a3(idim) * wmean(ixo^s,w_)**(adv_pow-1))
          cmax(ixo^s,1) = max(cmax(ixo^s,1), a3(idim) * wmean(ixo^s,w_)**(adv_pow-1))
       end if
    else
       cmax(ixo^s,1) = maxval(abs(a1(idim) * wmean(ixo^s,u_)**(adv_pow-1)))
       if (number_of_species >=2) then
          cmax(ixo^s,1) = max(cmax(ixo^s,1), maxval(abs(a2(idim) * wmean(ixo^s,v_)**(adv_pow-1))))
       end if
       if (number_of_species >=3) then
          cmax(ixo^s,1) = max(cmax(ixo^s,1), maxval(abs(a3(idim) * wmean(ixo^s,w_)**(adv_pow-1))))
       end if
    end if
 
  end subroutine ard_get_cbounds
 
  subroutine ard_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: dx^d, x(ixi^s, 1:^nd)
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(inout) :: dtnew
    double precision                :: maxrate
    double precision                :: maxd
    double precision                :: maxa
 
    dtnew = bigdouble
 
    ! dt < dx^2 / (2 * ndim * diffusion_coeff)
    ! use dtdiffpar < 1 for explicit and > 1 for imex/split
    maxd = d1
    if (number_of_species >= 2) then
        maxd = max(maxd, d2)
    end if
    if (number_of_species >= 3) then
        maxd = max(maxd, d3)
    end if
    dtnew = min(dtnew, dtdiffpar * minval([ dx^d ])**2 / (2 * ndim * maxd))
 
    ! Estimate time step for reactions
    select case (equation_type)
    case (eq_gray_scott)
       maxrate = max(maxval(w(ixo^s, v_))**2 + gs_f, &
            maxval(w(ixo^s, v_) * w(ixo^s, u_)) - gs_f - gs_k)
    case (eq_schnakenberg)
       maxrate = max(maxval(abs(w(ixo^s, v_) * w(ixo^s, u_) - 1)), &
            maxval(w(ixo^s, u_))**2)
    case (eq_brusselator)
       maxrate = max( maxval(w(ixo^s, u_)*w(ixo^s, v_) - (br_b+1)), &
            maxval(w(ixo^s, u_)**2) )
    case (eq_ext_brusselator)
       maxrate = max( maxval(w(ixo^s, u_)*w(ixo^s, v_) - (br_b+1)) + br_c, &
            maxval(w(ixo^s, u_)**2) )
       maxrate = max(maxrate, br_d)
    case (eq_logistic)
       maxrate = lg_lambda*maxval(abs(1 - w(ixo^s, u_))) ! abs for safety, normally u < 1
    case (eq_analyt_hunds)
       maxrate = maxval(w(ixo^s, u_)*abs(1 - w(ixo^s, u_))) / d1
    case (eq_belousov_fn)
       maxrate = max(&
            maxval(abs(1.0d0 - w(ixo^s, w_) - w(ixo^s, u_))) / bzfn_epsilon, &
            maxval(bzfn_lambda + w(ixo^s, u_)) / bzfn_delta &
       )
    case (eq_lorenz)
       ! det(J) = sigma(b(r-1) + x*(x*+y*))
       maxrate = max(lor_sigma, 1.0d0, lor_b)
    case (eq_no_reac)
       ! No reaction term, so no influence on timestep
       maxrate = zero
    case default
       maxrate = one
       call mpistop("Unknown equation type")
    end select
 
    dtnew = min(dtnew, dtreacpar / maxrate)
 
  end subroutine ard_get_dt
 
  ! Add the flux from the advection term
  subroutine ard_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wc(ixi^s, 1:nw)
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
    double precision, intent(out)   :: f(ixi^s, nwflux)
 
    f(ixo^s, u_) = (a1(idim)/adv_pow) * w(ixo^s,u_)**adv_pow
    if (number_of_species >=2) then
       f(ixo^s, v_) = (a2(idim)/adv_pow) * w(ixo^s,v_)**adv_pow
    end if
    if (number_of_species >=3) then
       f(ixo^s, w_) = (a3(idim)/adv_pow) * w(ixo^s,w_)**adv_pow
    end if
 
  end subroutine ard_get_flux
 
  subroutine ard_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim,w, x)
 
    ! Add geometrical source terms
    ! There are no geometrical source terms 
 
    use mod_global_parameters
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: qdt, dtfactor, x(ixi^s, 1:^nd)
    double precision, intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s,1:nw),w(ixi^s, 1:nw)
 
  end subroutine ard_add_source_geom
 
  ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
  subroutine ard_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, dtfactor
    double precision, intent(in)    :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw), x(ixi^s, 1:ndim)
    double precision, intent(inout) :: w(ixi^s, 1:nw)
    double precision                :: lpl_u(ixo^s), lpl_v(ixo^s), lpl_w(ixo^s)
    logical, intent(in)             :: qsourcesplit
    logical, intent(inout)          :: active
 
    ! here we add the reaction terms (always) and the diffusion if no imex is used
    if (qsourcesplit .eqv. ard_source_split) then
       if (.not.use_imex_scheme) then
          call ard_laplacian(ixi^l, ixo^l, wct(ixi^s, u_), lpl_u)
          if (number_of_species >= 2) then
             call ard_laplacian(ixi^l, ixo^l, wct(ixi^s, v_), lpl_v)
          end if
          if (number_of_species >= 3) then
             call ard_laplacian(ixi^l, ixo^l, wct(ixi^s, w_), lpl_w)
          end if
       else
          ! for all IMEX scheme variants: only add the reactions
          lpl_u = 0.0d0
          lpl_v = 0.0d0
          lpl_w = 0.0d0
       end if
 
       select case (equation_type)
       case (eq_gray_scott)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u - &
               wct(ixo^s, u_) * wct(ixo^s, v_)**2 + &
               gs_f * (1 - wct(ixo^s, u_)))
          w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v + &
               wct(ixo^s, u_) * wct(ixo^s, v_)**2 - &
               (gs_f + gs_k) * wct(ixo^s, v_))
       case (eq_schnakenberg)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + sb_kappa * (sb_alpha - wct(ixo^s, u_) + &
               wct(ixo^s, u_)**2 * wct(ixo^s, v_)))
          w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
               + sb_kappa * (sb_beta - wct(ixo^s, u_)**2 * wct(ixo^s, v_)))
       case (eq_brusselator)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + br_a - (br_b + 1) * wct(ixo^s, u_) &
               + wct(ixo^s, u_)**2 * wct(ixo^s, v_))
          w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
               + br_b * wct(ixo^s, u_) - wct(ixo^s, u_)**2 * wct(ixo^s, v_))
       case (eq_ext_brusselator)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + br_a - (br_b + 1) * wct(ixo^s, u_) &
               + wct(ixo^s, u_)**2 * wct(ixo^s, v_) &
               - br_c * wct(ixo^s, u_) + br_d * w(ixo^s, w_))
          w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
               + br_b * wct(ixo^s, u_) &
               - wct(ixo^s, u_)**2 * wct(ixo^s, v_))
          w(ixo^s, w_) = w(ixo^s, w_) + qdt * (d3 * lpl_w &
               + br_c * wct(ixo^s, u_) - br_d * w(ixo^s, w_))
       case (eq_logistic)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + lg_lambda * w(ixo^s, u_) * (1 - w(ixo^s, u_)))
       case (eq_analyt_hunds)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + 1.0d0/d1 * w(ixo^s, u_)**2 * (1 - w(ixo^s, u_)))
       case (eq_belousov_fn)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + 1.0d0/bzfn_epsilon * (bzfn_lambda * wct(ixo^s, u_) &
               - wct(ixo^s, u_)*wct(ixo^s, w_) + wct(ixo^s, u_) &
               - wct(ixo^s, u_)**2))
          w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
               + wct(ixo^s, u_) - wct(ixo^s, v_))
          w(ixo^s, w_) = w(ixo^s, w_) + qdt * (d3 * lpl_w &
               + 1.0d0/bzfn_delta * (-bzfn_lambda * wct(ixo^s, w_) &
               - wct(ixo^s, u_)*wct(ixo^s, w_) + bzfn_mu * wct(ixo^s, v_)))
       case (eq_lorenz)
          ! xdot = sigma.(y-x)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
               + lor_sigma * (wct(ixo^s, v_) - wct(ixo^s, u_)))
          ! ydot = r.x - y - x.z
          w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
               + lor_r * wct(ixo^s, u_) - wct(ixo^s, v_) &
               - wct(ixo^s, u_)*wct(ixo^s, w_))
          ! zdot = x.y - b.z
          w(ixo^s, w_) = w(ixo^s, w_) + qdt * (d3 * lpl_w &
               + wct(ixo^s, u_)*wct(ixo^s, v_) - lor_b * wct(ixo^s, w_))
       case (eq_no_reac)
          w(ixo^s, u_) = w(ixo^s, u_) + qdt * d1 * lpl_u
       case default
          call mpistop("Unknown equation type")
       end select
 
       ! enforce getbc call after source addition
       active = .true.
    end if
 
  end subroutine ard_add_source
 
  !> Compute the Laplacian using a standard second order scheme. For now this
  !> method only works in slab geometries. Requires one ghost cell only.
  subroutine ard_laplacian(ixI^L,ixO^L,var,lpl)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: var(ixi^s)
    double precision, intent(out) :: lpl(ixo^s)
    integer                       :: idir, jxo^l, hxo^l
    double precision              :: h_inv2
 
    if (slab) then
       lpl(ixo^s) = 0.0d0
       do idir = 1, ndim
          hxo^l=ixo^l-kr(idir,^d);
          jxo^l=ixo^l+kr(idir,^d);
          h_inv2 = 1/dxlevel(idir)**2
          lpl(ixo^s) = lpl(ixo^s) + h_inv2 * &
               (var(jxo^s) - 2 * var(ixo^s) + var(hxo^s))
       end do
    else
       call mpistop("ard_laplacian not implemented in this geometry")
    end if
 
  end subroutine ard_laplacian
 
  subroutine put_laplacians_onegrid(ixI^L,ixO^L,w)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, 1:nw)
 
    double precision                :: lpl_u(ixo^s), lpl_v(ixo^s), lpl_w(ixo^s)
 
    call ard_laplacian(ixi^l, ixo^l, w(ixi^s, u_), lpl_u)
    if (number_of_species >= 2) then
       call ard_laplacian(ixi^l, ixo^l, w(ixi^s, v_), lpl_v)
    end if
    if (number_of_species >= 3) then
       call ard_laplacian(ixi^l, ixo^l, w(ixi^s, w_), lpl_w)
    end if
 
    w(ixo^s,u_) = d1*lpl_u
    if (number_of_species >= 2) then
       w(ixo^s,v_) = d2*lpl_v
    end if
    if (number_of_species >= 3) then
       w(ixo^s,w_) = d3*lpl_w
    end if
 
  end subroutine put_laplacians_onegrid
  
  !> inplace update of psa==>F_im(psa)
  subroutine ard_evaluate_implicit(qtC,psa)
    use mod_global_parameters
    type(state), target :: psa(max_blocks)
    double precision, intent(in) :: qtc
 
    integer :: iigrid, igrid, level
    integer :: ixo^l
 
    !ixO^L=ixG^LL^LSUB1;
    ixo^l=ixm^ll;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       call put_laplacians_onegrid(ixg^ll,ixo^l,psa(igrid)%w)
    end do
    !$OMP END PARALLEL DO
 
  end subroutine ard_evaluate_implicit
 
  !> Implicit solve of psa=psb+dtfactor*dt*F_im(psa)
  subroutine ard_implicit_update(dtfactor,qdt,qtC,psa,psb)
    use mod_global_parameters
    use mod_forest
 
    type(state), target :: psa(max_blocks)
    type(state), target :: psb(max_blocks)
    double precision, intent(in) :: qdt
    double precision, intent(in) :: qtc
    double precision, intent(in) :: dtfactor
 
    integer                      :: n
    double precision             :: res, max_residual, lambda
 
    integer                        :: iw_to,iw_from
    integer                        :: iigrid, igrid, id
    integer                        :: nc, lvl
    type(tree_node), pointer       :: pnode
    real(dp)                       :: fac
 
    ! Avoid setting a very restrictive limit to the residual when the time step
    ! is small (as the operator is ~ 1/(D * qdt))
    if (qdt < dtmin) then
        if(mype==0)then
            print *,'skipping implicit solve: dt too small!'
            print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
        endif
        return
    endif
    max_residual = 1d-7/qdt
 
    mg%operator_type = mg_helmholtz
    call mg_set_methods(mg)
 
    if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
 
    ! First handle the u variable ***************************************
    lambda           = 1/(dtfactor * qdt * d1)
    call helmholtz_set_lambda(lambda)
    mg%bc(:, mg_iphi) = ard_mg_bc(1, :)
 
    call mg_copy_to_tree(u_, mg_irhs, factor=-lambda, state_from=psb)
    call mg_copy_to_tree(u_, mg_iphi, state_from=psb)
 
    call mg_fas_fmg(mg, .true., max_res=res)
    do n = 1, 10
       call mg_fas_vcycle(mg, max_res=res)
       if (res < max_residual) exit
    end do
 
    call mg_copy_from_tree_gc(mg_iphi, u_, state_to=psa)
    ! Done with the u variable ***************************************
 
    ! Next handle the v variable ***************************************
    if (number_of_species >= 2) then
       lambda = 1/(dtfactor * qdt * d2)
       call helmholtz_set_lambda(lambda)
       mg%bc(:, mg_iphi) = ard_mg_bc(2, :)
 
       call mg_copy_to_tree(v_, mg_irhs, factor=-lambda, state_from=psb)
       call mg_copy_to_tree(v_, mg_iphi, state_from=psb)
 
       call mg_fas_fmg(mg, .true., max_res=res)
       do n = 1, 10
          call mg_fas_vcycle(mg, max_res=res)
          if (res < max_residual) exit
       end do
 
       call mg_copy_from_tree_gc(mg_iphi, v_, state_to=psa)
    end if
    ! Done with the v variable ***************************************
 
    ! Next handle the w variable ***************************************
    if (number_of_species >= 3) then
       lambda = 1/(dtfactor * qdt * d3)
       call helmholtz_set_lambda(lambda)
 
       call mg_copy_to_tree(w_, mg_irhs, factor=-lambda, state_from=psb)
       call mg_copy_to_tree(w_, mg_iphi, state_from=psb)
 
       call mg_fas_fmg(mg, .true., max_res=res)
       do n = 1, 10
          call mg_fas_vcycle(mg, max_res=res)
          if (res < max_residual) exit
       end do
 
       call mg_copy_from_tree_gc(mg_iphi, w_, state_to=psa)
    end if
    ! Done with the w variable ***************************************
 
  end subroutine ard_implicit_update
 
end module mod_ard_phys