MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_rd_phys.t
Go to the documentation of this file.
1 !> Reaction-diffusion module (physics routines)
2 !>
3 !> Multiple reaction-diffusion systems are included: the Gray-Scott model, the
4 !> Schnakenberg model, the Brusselator model, the diffusive logistic equation,
5 !> an analytical testcase from "Numerical solution of time-dependent advection-
6 !> diffusion-reaction equations" by Hundsdorfer & Verwer, the Oregonator model,
7 !> the extended Brusselator model, and the diffusive Lorenz system. See the
8 !> documentation of the reaction-diffusion module for more information.
9 !>
10 !> IMEX methods are also supported. The implicit system is solved by a
11 !> multigrid solver coupled into MPI-AMRVAC.
12 !>
15  use mod_comm_lib, only: mpistop
16 
17 
18  implicit none
19  private
20 
21  integer, protected, public :: u_ = 1
22  integer, protected, public :: v_ = 2 !< For 2 or more equations
23  integer, protected, public :: w_ = 3 !< For 3 or more equations
24 
25  !> Whether particles module is added
26  logical, public, protected :: rd_particles = .false.
27 
28  !> Parameter with which to multiply the reaction timestep restriction
29  double precision, public, protected :: dtreacpar = 0.5d0
30 
31  !> Name of the system to be solved
32  character(len=20), public, protected :: equation_name = "gray-scott"
33  integer :: number_of_species = 2
34  integer :: equation_type = 1
35  integer, parameter :: eq_gray_scott = 1
36  integer, parameter :: eq_schnakenberg = 2
37  integer, parameter :: eq_brusselator = 3
38  integer, parameter :: eq_logistic = 4
39  integer, parameter :: eq_analyt_hunds = 5
40  integer, parameter :: eq_belousov_fn = 6 ! Field-Noyes model, or Oregonator
41  integer, parameter :: eq_ext_brusselator = 7
42  integer, parameter :: eq_lorenz = 8
43 
44  !> Diffusion coefficient for first species (u)
45  double precision, public, protected :: d1 = 0.05d0
46  !> Diffusion coefficient for second species (v) (if applicable)
47  double precision, public, protected :: d2 = 1.0d0
48  !> Diffusion coefficient for third species (w) (if applicable)
49  double precision, public, protected :: d3 = 1.0d0
50 
51  !> Parameter for Schnakenberg model
52  double precision, public, protected :: sb_alpha = 0.1305d0
53  !> Parameter for Schnakenberg model
54  double precision, public, protected :: sb_beta = 0.7695d0
55  !> Parameter for Schnakenberg model
56  double precision, public, protected :: sb_kappa = 100.0d0
57 
58  !> Feed rate for Gray-Scott model
59  double precision, public, protected :: gs_f = 0.046d0
60  !> Kill rate for Gray-Scott model
61  double precision, public, protected :: gs_k = 0.063d0
62 
63  !> Parameter for Brusselator model
64  double precision, public, protected :: br_a = 4.5d0
65  !> Parameter for Brusselator model
66  double precision, public, protected :: br_b = 8.0d0
67  !> Parameter for extended Brusselator model
68  double precision, public, protected :: br_c = 1.0d0
69  !> Parameter for extended Brusselator model
70  double precision, public, protected :: br_d = 1.0d0
71 
72  !> Parameter for logistic model (Fisher / KPP equation)
73  double precision, public, protected :: lg_lambda = 1.0d0
74 
75  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
76  double precision, public, protected :: bzfn_epsilon = 1.0d0
77  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
78  double precision, public, protected :: bzfn_delta = 1.0d0
79  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
80  double precision, public, protected :: bzfn_lambda = 1.0d0
81  !> Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction
82  double precision, public, protected :: bzfn_mu = 1.0d0
83 
84  !> Parameter for Lorenz system (Rayleigh number)
85  double precision, public, protected :: lor_r = 28.0d0
86  !> Parameter for Lorenz system (Prandtl number)
87  double precision, public, protected :: lor_sigma = 10.0d0
88  !> Parameter for Lorenz system (aspect ratio of the convection rolls)
89  double precision, public, protected :: lor_b = 8.0d0 / 3.0d0
90 
91  !> Whether to handle the explicitly handled source term in split fashion
92  logical :: rd_source_split = .false.
93 
94  !> Boundary condition information for the multigrid method
95  type(mg_bc_t), public :: rd_mg_bc(3, mg_num_neighbors)
96 
97  ! Public methods
98  public :: rd_phys_init
99 
100 contains
101 
102  subroutine rd_params_read(files)
103  use mod_global_parameters, only: unitpar
104  character(len=*), intent(in) :: files(:)
105  integer :: n
106 
107  namelist /rd_list/ d1, d2, d3, sb_alpha, sb_beta, sb_kappa, gs_f, gs_k, &
110  equation_name, rd_particles, rd_source_split, dtreacpar
111 
112  do n = 1, size(files)
113  open(unitpar, file=trim(files(n)), status='old')
114  read(unitpar, rd_list, end=111)
115 111 close(unitpar)
116  end do
117 
118  select case (equation_name)
119  case ("gray-scott")
120  equation_type = eq_gray_scott
121  number_of_species = 2
122  case ("schnakenberg")
123  equation_type = eq_schnakenberg
124  number_of_species = 2
125  case ("brusselator")
126  equation_type = eq_brusselator
127  number_of_species = 2
128  case ("logistic")
129  equation_type = eq_logistic
130  number_of_species = 1
131  case ("analyt_hunds")
132  equation_type = eq_analyt_hunds
133  number_of_species = 1
134  case ("belousov_fieldnoyes")
135  equation_type = eq_belousov_fn
136  number_of_species = 3
137  case ("ext_brusselator")
138  equation_type = eq_ext_brusselator
139  number_of_species = 3
140  case ("lorenz")
141  equation_type = eq_lorenz
142  number_of_species = 3
143  case default
144  call mpistop("Unknown equation_name (valid: gray-scott, schnakenberg, ...)")
145  end select
146 
147  end subroutine rd_params_read
148 
149  !> Write this module's parameters to a snapshot
150  subroutine rd_write_info(fh)
152  integer, intent(in) :: fh
153  integer, parameter :: n_par = 0
154  integer, dimension(MPI_STATUS_SIZE) :: st
155  integer :: er
156  integer :: idim
157 
158  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
159  end subroutine rd_write_info
160 
161  subroutine rd_phys_init()
163  use mod_physics
165  use mod_particles, only: particles_init
166 
167  call rd_params_read(par_files)
168 
169  physics_type = "rd"
170  phys_energy = .false.
171  phys_req_diagonal = .false.
173 
174  allocate(start_indices(number_species),stop_indices(number_species))
175  ! set the index of the first flux variable for species 1
176  start_indices(1)=1
177  ! Use the first variable as a density
178  u_ = var_set_fluxvar("u", "u")
179  if (number_of_species >= 2) then
180  v_ = var_set_fluxvar("v", "v")
181  end if
182  if (number_of_species >= 3) then
183  w_ = var_set_fluxvar("w", "w")
184  end if
185 
186  ! set number of variables which need update ghostcells
187  nwgc=nwflux
188 
189  ! set the index of the last flux variable for species 1
190  stop_indices(1)=nwflux
191 
192  ! Disable flux conservation near AMR boundaries, since we have no fluxes
193  fix_conserve_global = .false.
194 
206 
207  ! Initialize particles module
208  if (rd_particles) then
209  call particles_init()
210  phys_req_diagonal = .true.
211  end if
212 
213  end subroutine rd_phys_init
214 
215  subroutine rd_check_params
217  integer :: n, i, iw, species_list(number_of_species)
218 
219  if (any(flux_method /= fs_source)) then
220  ! there are no fluxes, only source terms in reaction-diffusion
221  call mpistop("mod_rd requires flux_scheme = source")
222  end if
223 
224  if (use_imex_scheme) then
225  use_multigrid = .true.
226  select case(number_of_species)
227  case(1); species_list = [u_]
228  case(2); species_list = [u_, v_]
229  case(3); species_list = [u_, v_, w_]
230  end select
231 
232  do i = 1, number_of_species
233  iw = species_list(i)
234 
235  ! Set boundary conditions for the multigrid solver
236  do n = 1, 2*ndim
237  select case (typeboundary(iw, n))
238  case (bc_symm)
239  ! d/dx u = 0
240  rd_mg_bc(i, n)%bc_type = mg_bc_neumann
241  rd_mg_bc(i, n)%bc_value = 0.0_dp
242  case (bc_asymm)
243  ! u = 0
244  rd_mg_bc(i, n)%bc_type = mg_bc_dirichlet
245  rd_mg_bc(i, n)%bc_value = 0.0_dp
246  case (bc_cont)
247  ! d/dx u = 0
248  rd_mg_bc(i, n)%bc_type = mg_bc_neumann
249  rd_mg_bc(i, n)%bc_value = 0.0_dp
250  case (bc_periodic)
251  ! Nothing to do here
252  case (bc_special)
253  if (.not. associated(rd_mg_bc(i, n)%boundary_cond)) then
254  write(*, "(A,I0,A,I0,A)") "typeboundary(", iw, ",", n, &
255  ") is 'special', but the corresponding method " // &
256  "rd_mg_bc(i, n)%boundary_cond is not set"
257  call mpistop("rd_mg_bc(i, n)%boundary_cond not set")
258  end if
259  case default
260  write(*,*) "rd_check_params warning: unknown boundary type"
261  rd_mg_bc(i, n)%bc_type = mg_bc_dirichlet
262  rd_mg_bc(i, n)%bc_value = 0.0_dp
263  end select
264  end do
265  end do
266  end if
267 
268  end subroutine rd_check_params
269 
270  subroutine rd_to_conserved(ixI^L, ixO^L, w, x)
272  integer, intent(in) :: ixI^L, ixO^L
273  double precision, intent(inout) :: w(ixI^S, nw)
274  double precision, intent(in) :: x(ixI^S, 1:^ND)
275 
276  ! Do nothing (primitive and conservative are equal for rd module)
277  end subroutine rd_to_conserved
278 
279  subroutine rd_to_primitive(ixI^L, ixO^L, w, x)
281  integer, intent(in) :: ixI^L, ixO^L
282  double precision, intent(inout) :: w(ixI^S, nw)
283  double precision, intent(in) :: x(ixI^S, 1:^ND)
284 
285  ! Do nothing (primitive and conservative are equal for rd module)
286  end subroutine rd_to_primitive
287 
288  subroutine rd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
290  integer, intent(in) :: ixI^L, ixO^L, idim
291  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S, 1:^ND)
292  double precision, intent(inout) :: cmax(ixI^S)
293 
294  cmax(ixo^s) = 0.0d0
295  end subroutine rd_get_cmax
296 
297  subroutine rd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed, cmax, cmin)
299  use mod_variables
300  integer, intent(in) :: ixI^L, ixO^L, idim
301  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S,nw)
302  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S,nw)
303  double precision, intent(in) :: x(ixI^S, 1:^ND)
304  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
305  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
306  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
307 
308  if (present(cmin)) then
309  cmin(ixo^s,1) = 0.0d0
310  cmax(ixo^s,1) = 0.0d0
311  else
312  cmax(ixo^s,1) = 0.0d0
313  end if
314 
315  end subroutine rd_get_cbounds
316 
317  subroutine rd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
319  integer, intent(in) :: ixI^L, ixO^L
320  double precision, intent(in) :: dx^D, x(ixI^S, 1:^ND)
321  double precision, intent(in) :: w(ixI^S, 1:nw)
322  double precision, intent(inout) :: dtnew
323  double precision :: maxrate
324  double precision :: maxD
325 
326  ! dt < dx^2 / (2 * ndim * diffusion_coeff)
327  ! use dtdiffpar < 1 for explicit and > 1 for imex/split
328  maxd = d1
329  if (number_of_species >= 2) then
330  maxd = max(maxd, d2)
331  end if
332  if (number_of_species >= 3) then
333  maxd = max(maxd, d3)
334  end if
335  dtnew = dtdiffpar * minval([ dx^d ])**2 / (2 * ndim * maxd)
336 
337  ! Estimate time step for reactions
338  select case (equation_type)
339  case (eq_gray_scott)
340  maxrate = max(maxval(w(ixo^s, v_))**2 + gs_f, &
341  maxval(w(ixo^s, v_) * w(ixo^s, u_)) - gs_f - gs_k)
342  case (eq_schnakenberg)
343  maxrate = max(maxval(abs(w(ixo^s, v_) * w(ixo^s, u_) - 1)), &
344  maxval(w(ixo^s, u_))**2)
345  case (eq_brusselator)
346  maxrate = max( maxval(w(ixo^s, u_)*w(ixo^s, v_) - (br_b+1)), &
347  maxval(w(ixo^s, u_)**2) )
348  case (eq_ext_brusselator)
349  maxrate = max( maxval(w(ixo^s, u_)*w(ixo^s, v_) - (br_b+1)) + br_c, &
350  maxval(w(ixo^s, u_)**2) )
351  maxrate = max(maxrate, br_d)
352  case (eq_logistic)
353  maxrate = lg_lambda*maxval(abs(1 - w(ixo^s, u_))) ! abs for safety, normally u < 1
354  case (eq_analyt_hunds)
355  maxrate = maxval(w(ixo^s, u_)*abs(1 - w(ixo^s, u_))) / d1
356  case (eq_belousov_fn)
357  maxrate = max(&
358  maxval(abs(1.0d0 - w(ixo^s, w_) - w(ixo^s, u_))) / bzfn_epsilon, &
359  maxval(bzfn_lambda + w(ixo^s, u_)) / bzfn_delta &
360  )
361  case (eq_lorenz)
362  ! det(J) = sigma(b(r-1) + x*(x*+y*))
363  maxrate = max(lor_sigma, 1.0d0, lor_b)
364  case default
365  maxrate = one
366  call mpistop("Unknown equation type")
367  end select
368 
369  dtnew = min(dtnew, dtreacpar / maxrate)
370 
371  end subroutine rd_get_dt
372 
373  ! There is nothing to add to the transport flux in the transport equation
374  subroutine rd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
376  integer, intent(in) :: ixI^L, ixO^L, idim
377  double precision, intent(in) :: wC(ixI^S, 1:nw)
378  double precision, intent(in) :: w(ixI^S, 1:nw)
379  double precision, intent(in) :: x(ixI^S, 1:^ND)
380  double precision, intent(out) :: f(ixI^S, nwflux)
381 
382  f(ixo^s, :) = 0.0d0
383  end subroutine rd_get_flux
384 
385  ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
386  subroutine rd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
388  integer, intent(in) :: ixI^L, ixO^L
389  double precision, intent(in) :: qdt, dtfactor
390  double precision, intent(in) :: wCT(ixI^S, 1:nw),wCTprim(ixI^S,1:nw),x(ixI^S, 1:ndim)
391  double precision, intent(inout) :: w(ixI^S, 1:nw)
392  double precision :: lpl_u(ixO^S), lpl_v(ixO^S), lpl_w(ixO^S)
393  logical, intent(in) :: qsourcesplit
394  logical, intent(inout) :: active
395 
396  ! here we add the reaction terms (always) and the diffusion if no imex is used
397  if (qsourcesplit .eqv. rd_source_split) then
398  if (.not.use_imex_scheme) then
399  call rd_laplacian(ixi^l, ixo^l, wct(ixi^s, u_), lpl_u)
400  if (number_of_species >= 2) then
401  call rd_laplacian(ixi^l, ixo^l, wct(ixi^s, v_), lpl_v)
402  end if
403  if (number_of_species >= 3) then
404  call rd_laplacian(ixi^l, ixo^l, wct(ixi^s, w_), lpl_w)
405  end if
406  else
407  ! for all IMEX scheme variants: only add the reactions
408  lpl_u = 0.0d0
409  lpl_v = 0.0d0
410  lpl_w = 0.0d0
411  end if
412 
413  select case (equation_type)
414  case (eq_gray_scott)
415  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u - &
416  wct(ixo^s, u_) * wct(ixo^s, v_)**2 + &
417  gs_f * (1 - wct(ixo^s, u_)))
418  w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v + &
419  wct(ixo^s, u_) * wct(ixo^s, v_)**2 - &
420  (gs_f + gs_k) * wct(ixo^s, v_))
421  case (eq_schnakenberg)
422  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
423  + sb_kappa * (sb_alpha - wct(ixo^s, u_) + &
424  wct(ixo^s, u_)**2 * wct(ixo^s, v_)))
425  w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
426  + sb_kappa * (sb_beta - wct(ixo^s, u_)**2 * wct(ixo^s, v_)))
427  case (eq_brusselator)
428  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
429  + br_a - (br_b + 1) * wct(ixo^s, u_) &
430  + wct(ixo^s, u_)**2 * wct(ixo^s, v_))
431  w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
432  + br_b * wct(ixo^s, u_) - wct(ixo^s, u_)**2 * wct(ixo^s, v_))
433  case (eq_ext_brusselator)
434  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
435  + br_a - (br_b + 1) * wct(ixo^s, u_) &
436  + wct(ixo^s, u_)**2 * wct(ixo^s, v_) &
437  - br_c * wct(ixo^s, u_) + br_d * w(ixo^s, w_))
438  w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
439  + br_b * wct(ixo^s, u_) &
440  - wct(ixo^s, u_)**2 * wct(ixo^s, v_))
441  w(ixo^s, w_) = w(ixo^s, w_) + qdt * (d3 * lpl_w &
442  + br_c * wct(ixo^s, u_) - br_d * w(ixo^s, w_))
443  case (eq_logistic)
444  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
445  + lg_lambda * w(ixo^s, u_) * (1 - w(ixo^s, u_)))
446  case (eq_analyt_hunds)
447  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
448  + 1.0d0/d1 * w(ixo^s, u_)**2 * (1 - w(ixo^s, u_)))
449  case (eq_belousov_fn)
450  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
451  + 1.0d0/bzfn_epsilon * (bzfn_lambda * wct(ixo^s, u_) &
452  - wct(ixo^s, u_)*wct(ixo^s, w_) + wct(ixo^s, u_) &
453  - wct(ixo^s, u_)**2))
454  w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
455  + wct(ixo^s, u_) - wct(ixo^s, v_))
456  w(ixo^s, w_) = w(ixo^s, w_) + qdt * (d3 * lpl_w &
457  + 1.0d0/bzfn_delta * (-bzfn_lambda * wct(ixo^s, w_) &
458  - wct(ixo^s, u_)*wct(ixo^s, w_) + bzfn_mu * wct(ixo^s, v_)))
459  case (eq_lorenz)
460  ! xdot = sigma.(y-x)
461  w(ixo^s, u_) = w(ixo^s, u_) + qdt * (d1 * lpl_u &
462  + lor_sigma * (wct(ixo^s, v_) - wct(ixo^s, u_)))
463  ! ydot = r.x - y - x.z
464  w(ixo^s, v_) = w(ixo^s, v_) + qdt * (d2 * lpl_v &
465  + lor_r * wct(ixo^s, u_) - wct(ixo^s, v_) &
466  - wct(ixo^s, u_)*wct(ixo^s, w_))
467  ! zdot = x.y - b.z
468  w(ixo^s, w_) = w(ixo^s, w_) + qdt * (d3 * lpl_w &
469  + wct(ixo^s, u_)*wct(ixo^s, v_) - lor_b * wct(ixo^s, w_))
470  case default
471  call mpistop("Unknown equation type")
472  end select
473 
474  ! enforce getbc call after source addition
475  active = .true.
476  end if
477 
478  end subroutine rd_add_source
479 
480  !> Compute the Laplacian using a standard second order scheme. For now this
481  !> method only works in slab geometries. Requires one ghost cell only.
482  subroutine rd_laplacian(ixI^L,ixO^L,var,lpl)
484  integer, intent(in) :: ixI^L, ixO^L
485  double precision, intent(in) :: var(ixI^S)
486  double precision, intent(out) :: lpl(ixO^S)
487  integer :: idir, jxO^L, hxO^L
488  double precision :: h_inv2
489 
490  if (slab) then
491  lpl(ixo^s) = 0.0d0
492  do idir = 1, ndim
493  hxo^l=ixo^l-kr(idir,^d);
494  jxo^l=ixo^l+kr(idir,^d);
495  h_inv2 = 1/dxlevel(idir)**2
496  lpl(ixo^s) = lpl(ixo^s) + h_inv2 * &
497  (var(jxo^s) - 2 * var(ixo^s) + var(hxo^s))
498  end do
499  else
500  call mpistop("rd_laplacian not implemented in this geometry")
501  end if
502 
503  end subroutine rd_laplacian
504 
505  subroutine put_laplacians_onegrid(ixI^L,ixO^L,w)
507  integer, intent(in) :: ixI^L, ixO^L
508  double precision, intent(inout) :: w(ixI^S, 1:nw)
509 
510  double precision :: lpl_u(ixO^S), lpl_v(ixO^S), lpl_w(ixO^S)
511 
512  call rd_laplacian(ixi^l, ixo^l, w(ixi^s, u_), lpl_u)
513  if (number_of_species >= 2) then
514  call rd_laplacian(ixi^l, ixo^l, w(ixi^s, v_), lpl_v)
515  end if
516  if (number_of_species >= 3) then
517  call rd_laplacian(ixi^l, ixo^l, w(ixi^s, w_), lpl_w)
518  end if
519 
520  w(ixo^s,u_) = d1*lpl_u
521  if (number_of_species >= 2) then
522  w(ixo^s,v_) = d2*lpl_v
523  end if
524  if (number_of_species >= 3) then
525  w(ixo^s,w_) = d3*lpl_w
526  end if
527 
528  end subroutine put_laplacians_onegrid
529 
530  !> inplace update of psa==>F_im(psa)
531  subroutine rd_evaluate_implicit(qtC,psa)
533  type(state), target :: psa(max_blocks)
534  double precision, intent(in) :: qtC
535 
536  integer :: iigrid, igrid, level
537  integer :: ixO^L
538 
539  !ixO^L=ixG^LL^LSUB1;
540  ixo^l=ixm^ll;
541  !$OMP PARALLEL DO PRIVATE(igrid)
542  do iigrid=1,igridstail; igrid=igrids(iigrid);
543  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
544  call put_laplacians_onegrid(ixg^ll,ixo^l,psa(igrid)%w)
545  end do
546  !$OMP END PARALLEL DO
547 
548  end subroutine rd_evaluate_implicit
549 
550  !> Implicit solve of psa=psb+dtfactor*dt*F_im(psa)
551  subroutine rd_implicit_update(dtfactor,qdt,qtC,psa,psb)
553  use mod_forest
554 
555  type(state), target :: psa(max_blocks)
556  type(state), target :: psb(max_blocks)
557  double precision, intent(in) :: qdt
558  double precision, intent(in) :: qtC
559  double precision, intent(in) :: dtfactor
560 
561  integer :: n
562  double precision :: res, max_residual, lambda
563 
564  integer :: iw_to,iw_from
565  integer :: iigrid, igrid, id
566  integer :: nc, lvl
567  type(tree_node), pointer :: pnode
568  real(dp) :: fac
569 
570  ! Avoid setting a very restrictive limit to the residual when the time step
571  ! is small (as the operator is ~ 1/(D * qdt))
572  if (qdt < dtmin) then
573  if(mype==0)then
574  print *,'skipping implicit solve: dt too small!'
575  print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
576  endif
577  return
578  endif
579  max_residual = 1d-7/qdt
580 
581  mg%operator_type = mg_helmholtz
582  call mg_set_methods(mg)
583 
584  if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
585 
586  ! First handle the u variable ***************************************
587  lambda = 1/(dtfactor * qdt * d1)
588  call helmholtz_set_lambda(lambda)
589  mg%bc(:, mg_iphi) = rd_mg_bc(1, :)
590 
591  call mg_copy_to_tree(u_, mg_irhs, factor=-lambda, state_from=psb)
592  call mg_copy_to_tree(u_, mg_iphi, state_from=psb)
593 
594  call mg_fas_fmg(mg, .true., max_res=res)
595  do n = 1, 10
596  call mg_fas_vcycle(mg, max_res=res)
597  if (res < max_residual) exit
598  end do
599 
600  call mg_copy_from_tree_gc(mg_iphi, u_, state_to=psa)
601  ! Done with the u variable ***************************************
602 
603  ! Next handle the v variable ***************************************
604  if (number_of_species >= 2) then
605  lambda = 1/(dtfactor * qdt * d2)
606  call helmholtz_set_lambda(lambda)
607  mg%bc(:, mg_iphi) = rd_mg_bc(2, :)
608 
609  call mg_copy_to_tree(v_, mg_irhs, factor=-lambda, state_from=psb)
610  call mg_copy_to_tree(v_, mg_iphi, state_from=psb)
611 
612  call mg_fas_fmg(mg, .true., max_res=res)
613  do n = 1, 10
614  call mg_fas_vcycle(mg, max_res=res)
615  if (res < max_residual) exit
616  end do
617 
618  call mg_copy_from_tree_gc(mg_iphi, v_, state_to=psa)
619  end if
620  ! Done with the v variable ***************************************
621 
622  ! Next handle the w variable ***************************************
623  if (number_of_species >= 3) then
624  lambda = 1/(dtfactor * qdt * d3)
625  call helmholtz_set_lambda(lambda)
626 
627  call mg_copy_to_tree(w_, mg_irhs, factor=-lambda, state_from=psb)
628  call mg_copy_to_tree(w_, mg_iphi, state_from=psb)
629 
630  call mg_fas_fmg(mg, .true., max_res=res)
631  do n = 1, 10
632  call mg_fas_vcycle(mg, max_res=res)
633  if (res < max_residual) exit
634  end do
635 
636  call mg_copy_from_tree_gc(mg_iphi, w_, state_to=psa)
637  end if
638  ! Done with the w variable ***************************************
639 
640  end subroutine rd_implicit_update
641 
642 end module mod_rd_phys
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module with basic grid data structures.
Definition: mod_forest.t:2
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer, dimension(:), allocatable flux_method
Which flux scheme of spatial discretization to use (per grid level)
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer, parameter bc_cont
integer, parameter bc_symm
logical fix_conserve_global
Whether to apply flux conservation at refinement boundaries.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
integer, parameter fs_source
double precision, dimension(ndim) dxlevel
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:15
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:59
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:81
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:68
procedure(sub_evaluate_implicit), pointer phys_evaluate_implicit
Definition: mod_physics.t:87
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:24
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:67
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:71
procedure(sub_check_params), pointer phys_check_params
Definition: mod_physics.t:56
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:58
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:16
procedure(sub_implicit_update), pointer phys_implicit_update
Definition: mod_physics.t:86
procedure(sub_add_source), pointer phys_add_source
Definition: mod_physics.t:73
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:61
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:27
Reaction-diffusion module (physics routines)
Definition: mod_rd_phys.t:13
subroutine rd_write_info(fh)
Write this module's parameters to a snapshot.
Definition: mod_rd_phys.t:151
double precision, public, protected gs_f
Feed rate for Gray-Scott model.
Definition: mod_rd_phys.t:59
double precision, public, protected lor_b
Parameter for Lorenz system (aspect ratio of the convection rolls)
Definition: mod_rd_phys.t:89
integer, public, protected u_
Definition: mod_rd_phys.t:21
double precision, public, protected bzfn_lambda
Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction.
Definition: mod_rd_phys.t:80
type(mg_bc_t), dimension(3, mg_num_neighbors), public rd_mg_bc
Boundary condition information for the multigrid method.
Definition: mod_rd_phys.t:95
subroutine, public rd_phys_init()
Definition: mod_rd_phys.t:162
subroutine rd_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Definition: mod_rd_phys.t:289
double precision, public, protected d3
Diffusion coefficient for third species (w) (if applicable)
Definition: mod_rd_phys.t:49
logical, public, protected rd_particles
Whether particles module is added.
Definition: mod_rd_phys.t:26
double precision, public, protected sb_kappa
Parameter for Schnakenberg model.
Definition: mod_rd_phys.t:56
double precision, public, protected br_d
Parameter for extended Brusselator model.
Definition: mod_rd_phys.t:70
double precision, public, protected gs_k
Kill rate for Gray-Scott model.
Definition: mod_rd_phys.t:61
double precision, public, protected d2
Diffusion coefficient for second species (v) (if applicable)
Definition: mod_rd_phys.t:47
subroutine rd_evaluate_implicit(qtC, psa)
inplace update of psa==>F_im(psa)
Definition: mod_rd_phys.t:532
double precision, public, protected dtreacpar
Parameter with which to multiply the reaction timestep restriction.
Definition: mod_rd_phys.t:29
integer, public, protected w_
For 3 or more equations.
Definition: mod_rd_phys.t:23
subroutine rd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Definition: mod_rd_phys.t:298
double precision, public, protected bzfn_epsilon
Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction.
Definition: mod_rd_phys.t:76
double precision, public, protected d1
Diffusion coefficient for first species (u)
Definition: mod_rd_phys.t:45
double precision, public, protected lg_lambda
Parameter for logistic model (Fisher / KPP equation)
Definition: mod_rd_phys.t:73
integer, public, protected v_
For 2 or more equations.
Definition: mod_rd_phys.t:22
double precision, public, protected sb_beta
Parameter for Schnakenberg model.
Definition: mod_rd_phys.t:54
subroutine rd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_rd_phys.t:318
double precision, public, protected bzfn_mu
Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction.
Definition: mod_rd_phys.t:82
double precision, public, protected sb_alpha
Parameter for Schnakenberg model.
Definition: mod_rd_phys.t:52
subroutine rd_check_params
Definition: mod_rd_phys.t:216
subroutine rd_to_primitive(ixIL, ixOL, w, x)
Definition: mod_rd_phys.t:280
subroutine rd_implicit_update(dtfactor, qdt, qtC, psa, psb)
Implicit solve of psa=psb+dtfactor*dt*F_im(psa)
Definition: mod_rd_phys.t:552
double precision, public, protected bzfn_delta
Parameter for the Field-Noyes model of the Belousov-Zhabotinsky reaction.
Definition: mod_rd_phys.t:78
double precision, public, protected lor_r
Parameter for Lorenz system (Rayleigh number)
Definition: mod_rd_phys.t:85
subroutine rd_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
Definition: mod_rd_phys.t:387
subroutine put_laplacians_onegrid(ixIL, ixOL, w)
Definition: mod_rd_phys.t:506
subroutine rd_to_conserved(ixIL, ixOL, w, x)
Definition: mod_rd_phys.t:271
double precision, public, protected br_a
Parameter for Brusselator model.
Definition: mod_rd_phys.t:64
subroutine rd_laplacian(ixIL, ixOL, var, lpl)
Compute the Laplacian using a standard second order scheme. For now this method only works in slab ge...
Definition: mod_rd_phys.t:483
double precision, public, protected lor_sigma
Parameter for Lorenz system (Prandtl number)
Definition: mod_rd_phys.t:87
subroutine rd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Definition: mod_rd_phys.t:375
double precision, public, protected br_c
Parameter for extended Brusselator model.
Definition: mod_rd_phys.t:68
character(len=20), public, protected equation_name
Name of the system to be solved.
Definition: mod_rd_phys.t:32
double precision, public, protected br_b
Parameter for Brusselator model.
Definition: mod_rd_phys.t:66
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11