Merge pull request #165 from AMReX-FHD/reactDiff

ReactDiff

exec/reactDiff/GNUmakefile

@@ -0,0 +1,64 @@
+# AMREX_HOME defines the directory in which we will find all the AMReX code.
+# If you set AMREX_HOME as an environment variable, this line will be ignored
+AMREX_HOME ?= ../../../amrex/
+
+DEBUG     = FALSE
+USE_MPI   = TRUE
+USE_OMP   = FALSE
+USE_CUDA  = FALSE
+COMP      = gnu
+DIM       = 2
+MAX_SPEC  = 8
+MAX_REAC  = 7
+
+TINY_PROFILE = FALSE
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+# add this back on if we add any local files
+#include ./Make.package
+#VPATH_LOCATIONS   += .
+#INCLUDE_LOCATIONS += .
+
+include ../../src_reactDiff/Make.package
+VPATH_LOCATIONS   += ../../src_reactDiff/
+INCLUDE_LOCATIONS += ../../src_reactDiff/
+
+include ../../src_chemistry/Make.package
+VPATH_LOCATIONS   += ../../src_chemistry/
+INCLUDE_LOCATIONS += ../../src_chemistry/
+
+include ../../src_analysis/Make.package
+VPATH_LOCATIONS   += ../../src_analysis/
+INCLUDE_LOCATIONS += ../../src_analysis/
+
+include ../../src_rng/Make.package
+VPATH_LOCATIONS   += ../../src_rng/
+INCLUDE_LOCATIONS += ../../src_rng/
+
+include ../../src_common/Make.package
+VPATH_LOCATIONS   += ../../src_common/
+INCLUDE_LOCATIONS += ../../src_common/
+
+include $(AMREX_HOME)/Src/Base/Make.package
+include $(AMREX_HOME)/Src/Boundary/Make.package
+include $(AMREX_HOME)/Src/LinearSolvers/MLMG/Make.package
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
+
+ifeq ($(findstring cgpu, $(HOST)), cgpu)
+  CXXFLAGS += $(FFTW)
+endif
+
+ifeq ($(USE_CUDA),TRUE)
+  LIBRARIES += -lcufft
+else
+  LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3
+endif
+
+MAXSPECIES := $(strip $(MAX_SPEC))
+DEFINES += -DMAX_SPECIES=$(MAXSPECIES)
+
+MAXREACTION := $(strip $(MAX_REAC))
+DEFINES += -DMAX_REACTION=$(MAXREACTION)
+

exec/reactDiff/inputs_paper_BPM_2d

@@ -0,0 +1,121 @@
+# This inputs file is used for generating
+# - Figures 5, 6, 7 (Section V.B)
+# in Paper by C. Kim et al. "Stochastic simulation of reaction-diffusion
+#  systems: A fluctuating-hydrodynamics approach"
+#  J. Chem. Phys. 146, 124110 (2017).
+# You can change some relevant parameters such as
+# - cell_depth
+# - n_cells (64^2 or 256^2) and max_grid_size
+# - fixed_dt, max_step
+# - plot_int (plot files)
+# - temporal_integrator, reaction_type (numerical scheme)
+# and run this inputs file.
+
+# Problem specification
+prob_lo =  0.0  0.0       # physical lo coordinate
+prob_hi = 32.0 32.0       # physical hi coordinate
+
+# number of cells in domain and maximum number of cells in a box
+n_cells = 64 64
+max_grid_size = 32 32
+
+# to compute cell volume in 2D problems
+cell_depth = 10.
+
+# Time-step control
+fixed_dt = 0.5
+
+# Controls for number of steps between actions
+max_step = 20000
+plot_int = 200
+struct_fact_int = -1
+n_steps_skip = 2000
+
+seed = 1
+
+nspecies = 3
+nreaction = 7
+
+prob_type = 0
+
+n_init_in_1 = 1685.8 533.5 56.38 # Start on the limit cycle
+
+integer_populations = 1
+
+# 0=D+R (first-order splitting)
+# 1=(1/2)R + D + (1/2)R (Strang option 1)
+# 2=(1/2)D + R + (1/2)D (Strang option 2)
+# -1=unsplit forward Euler
+# -2=unsplit explicit midpoint 
+# -3=unsplit multinomial diffusion
+# -4=unsplit implicit midpoint
+temporal_integrator = -4
+
+# only used for split schemes (temporal_integrator>=0)
+# 0=explicit trapezoidal predictor/corrector
+# 1=Crank-Nicolson semi-implicit
+# 2=explicit midpoint
+# 3=multinomial diffusion
+# 4=forward Euler  
+reactDiff_diffusion_type = 4
+
+# Fickian diffusion coeffs
+D_Fick = 0.1 0.01 0.01
+
+variance_coef_mass = 1.
+initial_variance_mass = 1.
+
+# how to compute n on faces for stochastic weighting
+# 1=arithmetic (with C0-Heaviside), 2=geometric, 3=harmonic
+# 10=arithmetic average with discontinuous Heaviside function
+# 11=arithmetic average with C1-smoothed Heaviside function
+# 12=arithmetic average with C2-smoothed Heaviside function
+avg_type = 1
+
+# only used for split schemes (temporal_integrator>=0)
+# 0=first-order (deterministic, tau leaping, CLE, or SSA)
+# 1=second-order (determinisitc, tau leaping, or CLE only)
+reactDiff_reaction_type = 0
+
+# 0=deterministic; 1=CLE; 2=SSA; 3=tau leap
+reaction_type = 3
+
+# BPM model is:
+#     (1) U + W --> V + W
+#     (2) V + V --> W
+#     (3) W     --> V + V
+#     (4) V     --> 0
+#     (5) 0     --> V
+#     (6) U     --> 0
+#     (7) 0     --> U
+stoich_1R = 1 0 1
+stoich_1P = 0 1 1
+stoich_2R = 0 2 0
+stoich_2P = 0 0 1
+stoich_3R = 0 0 1
+stoich_3P = 0 2 0
+stoich_4R = 0 1 0
+stoich_4P = 0 0 0
+stoich_5R = 0 0 0
+stoich_5P = 0 1 0
+stoich_6R = 1 0 0
+stoich_6P = 0 0 0
+stoich_7R = 0 0 0
+stoich_7P = 1 0 0
+
+# reaction rate constant for each reaction (assuming Law of Mass Action holds)
+# using rate_multiplier, reaction rates can be changed by the same factor
+# if include_discrete_LMA_correction, n^2 and n^3 in rate expressions become
+# n*(n-1/dv) and n*(n-1/dv)*(n-2/dv). 
+rate_const = 0.0002 0.0002 1. 0.03666663 4.44444555555 0.00333333 16.66665
+rate_multiplier = 1.
+include_discrete_LMA_correction = 1
+
+# Boundary conditions
+# ----------------------
+# BC specifications:
+# -1 = periodic
+#  1 = wall (Neumann)
+#  2 = reservoir (Dirichlet)
+bc_mass_lo = -1 -1
+bc_mass_hi = -1 -1

exec/reactDiff/inputs_paper_Lemarchand_3d

@@ -0,0 +1,109 @@
+# This inputs file is used for generating
+# - Figure 8 (Section V.C)
+# in Paper by C. Kim et al. "Stochastic simulation of reaction-diffusion
+#  systems: A fluctuating-hydrodynamics approach"
+#  J. Chem. Phys. 146, 124110 (2017).
+# You can change some relevant parameters such as
+# - initial_variance_mass: 0 (smooth initial condition) 1 (with fluctuations):
+# - variance_coef_mass: 0 (deterministic diffusion) 1 (stochastic)
+# - reaction_type: 0=deterministic; 1=CLE; 2=SSA; 3=tau leap
+# and run this inputs file.
+
+# Problem specification
+prob_lo =   0.0   0.0   0.0      # physical lo coordinate
+prob_hi = 512.0 512.0 512.0      # physical hi coordinate
+
+# number of cells in domain and maximum number of cells in a box
+n_cells = 256 256 256
+max_grid_size = 256 256 256
+
+# volume scale factor in 3D problems
+cell_depth = 1000.
+
+# Time-step control
+fixed_dt = 0.25
+
+# Controls for number of steps between actions
+max_step = 800
+plot_int = 10
+struct_fact_int = -1
+n_steps_skip = 200
+
+seed = 1
+
+nspecies = 2
+nreaction = 4
+
+prob_type = 5
+perturb_width = 16.    # scale factor for perturbed part in initial profile (for prob_type=4,5)
+smoothing_width = 1.   # scale factor for smoothing initial profile (for prob_type=4,5)
+
+n_init_in_1 = 2.16245 1.35018
+n_init_in_2 = 0. 10.
+
+# 0=D+R (first-order splitting)
+# 1=(1/2)R + D + (1/2)R (Strang option 1)
+# 2=(1/2)D + R + (1/2)D (Strang option 2)
+# -1=unsplit forward Euler
+# -2=unsplit explicit midpoint 
+# -3=unsplit multinomial diffusion
+# -4=unsplit implicit midpoint
+temporal_integrator = -4
+
+# only used for split schemes (temporal_integrator>=0)
+# 0=explicit trapezoidal predictor/corrector
+# 1=Crank-Nicolson semi-implicit
+# 2=explicit midpoint
+# 3=multinomial diffusion
+# 4=forward Euler  
+reactDiff_diffusion_type = 4
+
+# Fickian diffusion coeffs
+D_Fick = 1. 10.
+
+variance_coef_mass = 1.
+
+# how to compute n on faces for stochastic weighting
+# 1=arithmetic (with C0-Heaviside), 2=geometric, 3=harmonic
+# 10=arithmetic average with discontinuous Heaviside function
+# 11=arithmetic average with C1-smoothed Heaviside function
+# 12=arithmetic average with C2-smoothed Heaviside function
+avg_type = 1
+
+# only used for split schemes (temporal_integrator>=0)
+# 0=first-order (deterministic, tau leaping, CLE, or SSA)
+# 1=second-order (determinisitc, tau leaping, or CLE only)
+reactDiff_reaction_type = 0
+
+# 0=deterministic; 1=CLE; 2=SSA; 3=tau leap
+reaction_type = 3
+
+#     (1) A -> 0
+#     (2) 2A + B --> 3A
+#     (3) B --> 0
+#     (4) 0 --> B
+stoich_1R = 1 0
+stoich_1P = 0 0
+stoich_2R = 2 1
+stoich_2P = 3 0
+stoich_3R = 0 1
+stoich_3P = 0 0
+stoich_4R = 0 0
+stoich_4P = 0 1
+
+# reaction rate constant for each reaction (assuming Law of Mass Action holds)
+# using rate_multiplier, reaction rates can be changed by the same factor
+# if include_discrete_LMA_correction, n^2 and n^3 in rate expressions become
+# n*(n-1/dv) and n*(n-1/dv)*(n-2/dv). 
+rate_const =  4. 1.37 1. 10.
+rate_multiplier = 0.1
+include_discrete_LMA_correction = 1
+
+# Boundary conditions
+# ----------------------
+# BC specifications:
+# -1 = periodic
+#  1 = wall (Neumann)
+#  2 = reservoir (Dirichlet)
+bc_mass_lo = -1 -1 -1
+bc_mass_hi = -1 -1 -1

src_chemistry/chemistry_functions.H

@@ -19,4 +19,16 @@ void InitializeChemistryNamespace();
 // used in compressible code only
 void compute_compressible_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
                                                MultiFab& prim, MultiFab& source, MultiFab& ranchem);
+
+void ChemicalRates(const MultiFab& n_cc, MultiFab& chem_rate, const amrex::Geometry& geom, const amrex::Real& dt, 
+                   const MultiFab& n_interm, Vector<Real>& lin_comb_coef_in, Real volume_factor_in=1.);
+
+AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<Real,MAX_SPECIES>& n_in,
+                                                  GpuArray<Real,MAX_REACTION>& reaction_rates,
+                                                  const amrex::Real& dv);
+
+AMREX_GPU_HOST_DEVICE void sample_num_reactions(GpuArray<Real,MAX_SPECIES>& n_in,
+                                                GpuArray<Real,MAX_REACTION>& num_reactions,
+                                                GpuArray<Real,MAX_REACTION>& avg_num_reactions,
+                                                const amrex::RandomEngine& engine);
 #endif