mod_rd_phys.t

Go to the documentation of this file.

!> Reaction-diffusion module (physics routines)
!>
!> Multiple reaction-diffusion systems are included: 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, and the diffusive Lorenz system. See the
!> documentation of the reaction-diffusion module for more information.
!>
!> IMEX methods are also supported. The implicit system is solved by a
!> multigrid solver coupled into MPI-AMRVAC.
!>
module mod_rd_phys
  use mod_multigrid_coupling
  use mod_comm_lib, only: mpistop
 
 
  implicit none
  private
 
  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 :: rd_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
  integer, parameter :: eq_schnakenberg    = 2
  integer, parameter :: eq_brusselator     = 3
  integer, parameter :: eq_logistic        = 4
  integer, parameter :: eq_analyt_hunds    = 5
  integer, parameter :: eq_belousov_fn     = 6 ! Field-Noyes model, or Oregonator
  integer, parameter :: eq_ext_brusselator = 7
  integer, parameter :: eq_lorenz          = 8
 
  !> 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
 
  !> Parameter for Schnakenberg model
  double precision, public, protected :: sb_alpha = 0.1305d0
  !> Parameter for Schnakenberg model
  double precision, public, protected :: sb_beta  = 0.7695d0
  !> Parameter for Schnakenberg model
  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
 
  !> Parameter for Brusselator model
  double precision, public, protected :: br_a = 4.5d0
  !> Parameter for Brusselator model
  double precision, public, protected :: br_b = 8.0d0
  !> Parameter for extended Brusselator model
  double precision, public, protected :: br_c = 1.0d0
  !> Parameter for extended Brusselator model
  double precision, public, protected :: br_d = 1.0d0
 
  !> Parameter for logistic model (Fisher / KPP equation)
  double precision, public, protected :: lg_lambda = 1.0d0
 
  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
  double precision, public, protected :: bzfn_epsilon = 1.0d0
  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
  double precision, public, protected :: bzfn_delta   = 1.0d0
  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
  double precision, public, protected :: bzfn_lambda  = 1.0d0
  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
  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 :: rd_source_split = .false.
 
  !> Boundary condition information for the multigrid method
  type(mg_bc_t), public :: rd_mg_bc(3, mg_num_neighbors)
 
  ! Public methods
  public :: rd_phys_init
 
contains
 
  subroutine rd_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /rd_list/ d1, d2, d3, 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, rd_particles, rd_source_split, dtreacpar
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status='old')
       read(unitpar, rd_list, end=111)
111    close(unitpar)
    end do
 
    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 default
       call mpistop("Unknown equation_name (valid: gray-scott, schnakenberg, ...)")
    end select
 
  end subroutine rd_params_read
 
  !> Write this module's parameters to a snapshot
  subroutine rd_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 rd_write_info
 
  subroutine rd_phys_init()
    use mod_global_parameters
    use mod_physics
    use mod_multigrid_coupling
    use mod_particles, only: particles_init
 
    call rd_params_read(par_files)
 
    physics_type = "rd"
    phys_energy  = .false.
    use_particles = rd_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
  
    !  Number of variables need reconstruction in w
    nw_recon=nwflux
 
    ! Disable flux conservation near AMR boundaries, since we have no fluxes
    fix_conserve_global = .false.
 
    phys_get_cmax     => rd_get_cmax
    phys_get_cbounds  => rd_get_cbounds
    phys_get_flux     => rd_get_flux
    phys_to_conserved => rd_to_conserved
    phys_to_primitive => rd_to_primitive
    phys_add_source   => rd_add_source
    phys_get_dt       => rd_get_dt
    phys_write_info   => rd_write_info
    phys_check_params => rd_check_params
    phys_implicit_update   => rd_implicit_update
    phys_evaluate_implicit => rd_evaluate_implicit
 
    ! Initialize particles module
    if (rd_particles) then
       call particles_init()
    end if
 
  end subroutine rd_phys_init
 
  subroutine rd_check_params
    use mod_global_parameters
    integer :: n, i, iw, species_list(number_of_species)
 
    if (any(flux_method /= fs_source)) then
       ! there are no fluxes, only source terms in reaction-diffusion
       call mpistop("mod_rd requires flux_scheme = source")
    end if
 
    if (use_imex_scheme) then
       use_multigrid = .true.
       select case(number_of_species)
       case(1); species_list = [u_]
       case(2); species_list = [u_, v_]
       case(3); species_list = [u_, v_, w_]
       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
                rd_mg_bc(i, n)%bc_type = mg_bc_neumann
                rd_mg_bc(i, n)%bc_value = 0.0_dp
             case (bc_asymm)
                ! u = 0
                rd_mg_bc(i, n)%bc_type = mg_bc_dirichlet
                rd_mg_bc(i, n)%bc_value = 0.0_dp
             case (bc_cont)
                ! d/dx u = 0
                rd_mg_bc(i, n)%bc_type = mg_bc_neumann
                rd_mg_bc(i, n)%bc_value = 0.0_dp
             case (bc_periodic)
                ! Nothing to do here
             case (bc_special)
                if (.not. associated(rd_mg_bc(i, n)%boundary_cond)) then
                   write(*, "(A,I0,A,I0,A)") "typeboundary(", iw, ",", n, &
                        ") is 'special', but the corresponding method " // &
                        "rd_mg_bc(i, n)%boundary_cond is not set"
                   call mpistop("rd_mg_bc(i, n)%boundary_cond not set")
                end if
             case default
                write(*,*) "rd_check_params warning: unknown boundary type"
                rd_mg_bc(i, n)%bc_type = mg_bc_dirichlet
                rd_mg_bc(i, n)%bc_value = 0.0_dp
             end select
          end do
       end do
    end if
 
  end subroutine rd_check_params
 
  subroutine rd_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 rd module)
  end subroutine rd_to_conserved
 
  subroutine rd_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 rd module)
  end subroutine rd_to_primitive
 
  subroutine rd_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) = 0.0d0
  end subroutine rd_get_cmax
 
  subroutine rd_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(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
 
    if (present(cmin)) then
       cmin(ixo^s,1) = 0.0d0
       cmax(ixo^s,1) = 0.0d0
    else
       cmax(ixo^s,1) = 0.0d0
    end if
 
  end subroutine rd_get_cbounds
 
  subroutine rd_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
 
    ! 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 = 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 default
       maxrate = one
       call mpistop("Unknown equation type")
    end select
 
    dtnew = min(dtnew, dtreacpar / maxrate)
 
  end subroutine rd_get_dt
 
  ! There is nothing to add to the transport flux in the transport equation
  subroutine rd_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, :) = 0.0d0
  end subroutine rd_get_flux
 
  ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
  subroutine rd_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. rd_source_split) then
       if (.not.use_imex_scheme) then
          call rd_laplacian(ixi^l, ixo^l, wct(ixi^s, u_), lpl_u)
          if (number_of_species >= 2) then
             call rd_laplacian(ixi^l, ixo^l, wct(ixi^s, v_), lpl_v)
          end if
          if (number_of_species >= 3) then
             call rd_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 default
          call mpistop("Unknown equation type")
       end select
 
       ! enforce getbc call after source addition
       active = .true.
    end if
 
  end subroutine rd_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 rd_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("rd_laplacian not implemented in this geometry")
    end if
 
  end subroutine rd_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 rd_laplacian(ixi^l, ixo^l, w(ixi^s, u_), lpl_u)
    if (number_of_species >= 2) then
       call rd_laplacian(ixi^l, ixo^l, w(ixi^s, v_), lpl_v)
    end if
    if (number_of_species >= 3) then
       call rd_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 rd_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 rd_evaluate_implicit
 
  !> Implicit solve of psa=psb+dtfactor*dt*F_im(psa)
  subroutine rd_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) = rd_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) = rd_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 rd_implicit_update
 
end module mod_rd_phys