diff --git a/staggered_grid/exec/chargedFluid/test/Makefile.Donev b/staggered_grid/exec/chargedFluid/test/Makefile.Donev new file mode 100644 index 00000000..f7712e65 --- /dev/null +++ b/staggered_grid/exec/chargedFluid/test/Makefile.Donev @@ -0,0 +1,103 @@ +# FBOXLIB_HOME defines the directory in which we will find all the amrex code +# If you set FBOXLIB_HOME as an environment variable, this line will be ignored +FBOXLIB_HOME := ../../../../../FBoxLib +HYDROLIB_HOME := ../../../../HydroGrid +STAGGERED_HOME := ../../.. + +# Instructions: +# NDEBUG=t means optimize, not setting it means debug mode +# MPI=t means use MPI, otherwise compile for serial execution +# Important: Choose compilers (gfortran or ifort) + +NDEBUG := t +MPI := +OMP := +PROF := +COMP := gfortran +CCOMP := gcc +MKVERBOSE := t + +# For Intel compiler (requires ifx+icc as of version 2023) +#COMP := Intel +#CCOMP := icc + +# need this to compile bl_rng.f90 (module to store random number engines and distributions) +CXX11 = t + +include $(FBOXLIB_HOME)/Tools/F_mk/GMakedefs.mak + +include $(STAGGERED_HOME)/src_lowMach/GPackage.mak +VPATH_LOCATIONS += $(STAGGERED_HOME)/src_lowMach + +include $(STAGGERED_HOME)/src_charged/GPackage.mak +VPATH_LOCATIONS += $(STAGGERED_HOME)/src_charged + +include $(STAGGERED_HOME)/src_multiSpec/GPackage.mak +VPATH_LOCATIONS += $(STAGGERED_HOME)/src_multiSpec + +include $(STAGGERED_HOME)/src_chemistry/GPackage.mak +VPATH_LOCATIONS += $(STAGGERED_HOME)/src_chemistry + +include $(STAGGERED_HOME)/src_gmres/GPackage.mak +VPATH_LOCATIONS += $(STAGGERED_HOME)/src_gmres + +include $(STAGGERED_HOME)/src_common/GPackage.mak +VPATH_LOCATIONS += $(STAGGERED_HOME)/src_common + +include $(HYDROLIB_HOME)/GPackage.mak +VPATH_LOCATIONS += $(HYDROLIB_HOME) + +include $(FBOXLIB_HOME)/Src/BaseLib/GPackage.mak +VPATH_LOCATIONS += $(FBOXLIB_HOME)/Src/BaseLib + +include $(FBOXLIB_HOME)/Src/MultiGrid/GPackage.mak +VPATH_LOCATIONS += $(FBOXLIB_HOME)/Src/MultiGrid/ + +# Because module files are not portable, we would need to recompile +# LAPACK95 unless we use Intel for everything. +# Therefore, Donev has disabled LAPACK95 +#LAPACK95_DIR = $(STAGGERED_HOME)/src_multiSpec/LAPACK95 +#F90FLAGS += -I$(LAPACK95_DIR)/lapack95_modules + +# Linking lapack and blas, either default ones or Intel ones +# Intel OneAPI is assumed to be installed, which it is on both CIMS and donevm laptop +# For Intel, best is to just link this dynamic "use all of MKL" +# $(MKLROOT)/lib/intel64/libmkl_rt.so +# If you want static then link in +# $(MKLROOT)/lib/intel64/libmkl_blas95_lp64.a +# $(MKLROOT)/lib/intel64/libmkl_lapack95_lp64.a + +ifeq ($(findstring cims.nyu.edu, $(HOSTNAME)), cims.nyu.edu) + # Best to use Intel's LAPACK/BLAS on CIMS to avoid old gcc versions in the mix + # Make sure to "source /opt/pkg/intel/oneapi/setvars.sh" for ifort (module load intel-oneapi) + # This exports variables $(MKLROOT) and $(I_MPI_ROOT) + #MKLROOT=/opt/pkg/intel/oneapi/mkl/latest +else ifeq ($(findstring pop-os, $(HOSTNAME)), pop-os) + # ubuntu system76 pop-os + # Make sure to "source /home/donev/intel/oneapi/setvars.sh" for ifort + # MKLROOT=/home/donev/intel/oneapi/mkl/latest +endif + +ifeq ($(findstring gfortran, $(FC)), gfortran) + # Since LAPACK/BLAS are compiled with old gcc version + # Use instead the Intel MPI library + ifdef MKLROOT + # Use the single-library dynamic interface of Intel meant for GNU compiler: + libraries += $(MKLROOT)/lib/intel64/libmkl_rt.so + MKL_INTERFACE_LAYER=MKL_INTERFACE_LP64+MKL_INTERFACE_GNU + else + libraries += -llapack -lblas + endif +else ifeq ($(findstring if, $(FC)), if) # ifx or ifort (Intel) + # The Intel compiler mkl flag should do everything on its own, including linking: + F90FLAGS += -qmkl + MKL_INTERFACE_LAYER=MKL_INTERFACE_LP64 +else + libraries += -llapack -lblas +endif + +main.$(suf).exe: $(objects) + $(LINK.f90) -o main.$(suf).exe $(objects) $(libraries) + +include $(FBOXLIB_HOME)/Tools/F_mk/GMakerules.mak + diff --git a/staggered_grid/exec/chargedFluid/test/README b/staggered_grid/exec/chargedFluid/test/README index 573fcf0e..547e032d 100644 --- a/staggered_grid/exec/chargedFluid/test/README +++ b/staggered_grid/exec/chargedFluid/test/README @@ -1,3 +1,11 @@ +NEW -- two directories + +electro_osmosis_case1 +electro_osmosis_case2 + +are examples that run by reading in a pnp profile + + Supported Problems (./inputs_files/): ---------------------------------- diff --git a/staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/debye_huckel_check.py b/staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/debye_huckel_check.py new file mode 100644 index 00000000..0936e685 --- /dev/null +++ b/staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/debye_huckel_check.py @@ -0,0 +1,107 @@ +import numpy as np +import matplotlib.pyplot as plt + +debye_len_case = '1_50' +Ly = 1.28e-4 +#debye_len_case = '1_10' + +if (debye_len_case == '1_50'): + E_ext = 3.e9 + dphi_dn_wall = 794201762623.94604 +elif (debye_len_case == '1_10'): + E_ext = 3.e10 + dphi_dn_wall = 31767657979.612209 + +# import deterministic potential and velocity + +u_det = np.loadtxt('u_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +epot = np.loadtxt('epot_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +c1 = np.loadtxt('c1_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +c2 = np.loadtxt('c2_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +c3 = np.loadtxt('c3_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) + + +##################### +## need to calculuate lambda_d based on centerline concentration values: +##################### +mu = 1.05e-2 # momentum diffusion coefficient +dielectric_const = 6.91e-19 # dielectric constant +k_B = 1.3806488e-16 +T = 300. + +q1 = 4.2e3 +q2 = -2.72e3 + +m1 = 3.82e-23 +m2 = 5.89e-23 +m3 = 3.35e-23 +rho = 1. +c1_c = c1[len(c1)/2] +c2_c = c2[len(c2)/2] +c3_c = c3[len(c3)/2] +epot_c = epot[len(epot)/2] +u_c = u_det[len(u_det)/2] + +###################### +## check concentration profile against Boltzmann distribution +###################### +c1_boltz = c1_c*np.exp(-m1*q1/k_B/T*(epot-epot_c)) +c1_approx = c1_c*(1. - m1*q1/k_B/T*(epot-epot_c)) + +c2_boltz = c2_c*np.exp(-m2*q2/k_B/T*(epot-epot_c)) +c2_approx = c2_c*(1. - m2*q2/k_B/T*(epot-epot_c)) + + +##################### +## check relationship: +## +## u = \epsilon E_x/ \mu (\phi - \phi_{wall}) +## +##################### +shift = epot_c - u_c*mu/dielectric_const/E_ext +epot_wall = epot[0] +u_theory = dielectric_const*E_ext/mu * (epot - shift) + +###################### +## compute debye length based on concentrations at channel center +###################### +lambda_d = np.sqrt(dielectric_const*k_B*T/(rho*(c1_c*m1*q1**2 + c2_c*m2*q2**2))) +print 'debye len is: ', lambda_d +print 'ratio Ly/lambda_d: ', Ly/lambda_d + + + + +###################### +## compare centerline velocity with D-H theory +###################### +print 'Debye Huckel slip velocity: ', lambda_d*dielectric_const*dphi_dn_wall*E_ext/mu +print 'Deterministic calculation vel: ', u_c + + + +plt.figure(1) +plt.plot(u_det, 'b-o', label = 'simulation $u$') +plt.plot(u_theory, 'r-o', label = '$\epsilon E_x (\phi-\phi_w)/\mu$') +plt.title('Velocity') +plt.legend(loc = 'best') + +plt.figure(5) +plt.plot(u_det-u_theory, 'g-o') +plt.title('Difference between velocities') + +plt.figure(2) +plt.plot(c1, 'g-o', label = 'simulation') +plt.plot(c1_boltz, 'r-o', label = 'theory') +plt.plot(c1_approx, 'b-o', label = 'approximation') +plt.title('Concentration 1') +plt.legend(loc='best') + +plt.figure(3) +plt.plot(c2, 'g-o', label = 'simulation') +plt.plot(c2_boltz, 'r-o', label = 'theory') +plt.plot(c2_approx, 'b-o', label = 'approximation') +plt.title('Concentration 2') +plt.legend(loc='best') + +plt.show() diff --git a/staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/plot_1d_profs.py b/staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/plot_1d_profs.py new file mode 100644 index 00000000..09eeb086 --- /dev/null +++ b/staggered_grid/exec/chargedFluid/test/eo_post_proc_python_scripts/plot_1d_profs.py @@ -0,0 +1,76 @@ +import os +import urllib2 +import numpy as np +import matplotlib.font_manager as font_manager +import matplotlib.pyplot as plt +import pylab + + +##========================================== +## Get time averaged profiles +## + +#debye_len_case = '1_50' +Ly = 1.28e-4 +debye_len_case = '1_10' + +if (debye_len_case == '1_50'): + y_vals, u, v, w, c1, c2, c3, epot = np.loadtxt('hstat00064000_debye_len_1_50', unpack = True, usecols = [0,1,2,3,5,6,7,9]) +elif (debye_len_case == '1_10'): + y_vals, u, v, w, c1, c2, c3, epot = np.loadtxt('hstat00056500_debye_len_1_10', unpack = True, usecols = [0,1,2,3,5,6,7,9]) + +u_det = np.loadtxt('u_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +c1_det = np.loadtxt('c1_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +c2_det = np.loadtxt('c2_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +c3_det = np.loadtxt('c3_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) +epot_det = np.loadtxt('epot_deterministic_' + debye_len_case + '.txt', unpack = True, usecols = [0]) + + +##========================================== +## Plot results +## + +plt.figure(1) +plt.plot(y_vals, u, 'b-o', label = '3d averaged') +plt.plot(y_vals, u_det, 'r-o', label = '2d') +plt.title('$u$ v. wall-normal $y$') +plt.legend(loc='best') + +plt.figure(7) +plt.plot(y_vals, v, 'b-o', label = '3d averaged') +plt.title('$v$ v. wall-normal $y$') +plt.legend(loc='best') + +plt.figure(9) +plt.plot(y_vals, w, 'b-o', label = '3d averaged') +plt.title('$w$ v. wall-normal $y$') +plt.legend(loc='best') + +plt.figure(2) +plt.plot(y_vals, c1, 'b-o', label = '3d averaged') +plt.plot(y_vals, c1_det, 'r-o', label = '2d') +plt.title('$c_1$ v. wall-normal $y$') +plt.legend(loc='best') + +plt.figure(3) +plt.plot(y_vals, c2, 'b-o', label = '3d averaged') +plt.plot(y_vals, c2_det, 'r-o', label = '2d') +plt.title('$c_2$ v. wall-normal $y$') +plt.legend(loc='best') + +plt.figure(4) +plt.plot(y_vals, c3, 'b-o', label = '3d averaged') +plt.plot(y_vals, c3_det, 'r-o', label = '2d') +plt.title('$c_3$ v. wall-normal $y$') +plt.legend(loc='best') + +plt.figure(5) +plt.plot(y_vals, epot, 'b-o', label = '3d averaged') +plt.plot(y_vals, epot_det, 'r-o', label = '2d') +plt.title('$\phi$ v. wall-normal $y$') +plt.legend(loc='best') + +### show ### +plt.show() + + diff --git a/staggered_grid/exec/chargedFluid/test/inputs_arefh b/staggered_grid/exec/chargedFluid/test/inputs_arefh new file mode 100644 index 00000000..53f3c90a --- /dev/null +++ b/staggered_grid/exec/chargedFluid/test/inputs_arefh @@ -0,0 +1,303 @@ +! Options for low Mach / Boussinesq equations +&probin_common + + ! Geometry specification + !---------------------- + dim_in = 2 ! 2D or 3D + prob_lo(1:3) = 0.d0 0.d0 0.d0 ! physical lo coordinate + prob_hi(1:3) = 4.0 128.0 ! physical hi coordinate = [Nx*h, Ny*h] + ! Here dx=dy=h, and n_cells=[Nx, Ny] + n_cells(1:2) = 4 128 ! number of cells in domain + ! max number of cells in a box: + !max_grid_size(1:2) = 4 1024 ! Disable parallelization + max_grid_size(1:2) = 4 32 ! Parallelize into 4 boxes = 128/32 + + prob_type = 2 ! Sets up physical problem (see src_lowMach/init_lowmach.f90) + ! case 2: + ! constant concentration gradient along y + ! c=c_init(1,:) on bottom, c=c_init(2,:) on top + !---------- + ! case 17: For electro-osmotic and other channel flows (2D and 3D) + ! Concentrations and streamwise vel. u only depend on wall normal variable y, + ! ie the fields are homogeneous in the x and z directions. + ! The values for the fields are read in from .txt files. Also, v = w = 0. + ! Files are: "c1_1d_vals.txt", "c2_1d_vals.txt", "c3_1d_vals.txt", "u_1d_vals.txt" + ! Note that currently only nspecies=3 is supported for prob_type=17 + + + ! Grid size and time-step control + !---------------------- + fixed_dt = 0.05 ! time step (if positive, fixed) + max_step = 100 ! maximum number of time steps + plot_int = 1000 ! Interval for writing a plotfile (for visit/amrvis) + + ! Controls for number of steps between actions (for HydroGrid see below) + !---------------------- + chk_int = -1 ! Interval for writing a checkpoint + restart = -1 ! checkpoint restart number + print_int = 10 ! how often to output EOS drift and sum of conserved quantities + project_eos_int = 100 ! how often to call project_onto_eos + + ! Physical parameters + !-------------------- + nspecies = 3 ! Number of species (two ions and water) + molmass(1:3) = 1.0 1.0 1.0 ! molecular masses (in some units, make all equal for simplicity) + rhobar(1:3) = 1.d0 1.d0 1.d0 ! density, keep at 1, only used in low Mach codes, not here + grav(1:2) = 0.d0 0.d0 ! gravity vector (negative is downwards) + + ! stochastic forcing amplitudes (1 for physical values, 0 to run them off) + variance_coef_mom = 0.d0 ! Multiplier for stochastic momentum flux + variance_coef_mass = 0.d0 ! Multiplier for stochastic mass flux + k_B = 1 ! Boltzmann's constant, true = 1.38064852d-16 + + + ! Algorithm control / selection + !---------------------- + algorithm_type = 6 ! Select type of temporal integration + ! In low Mach codes: + ! 0 = Inertial trapezoidal lowMach + ! 2 = Overdamped lowMach with 2 RNGs + ! 3 = Iterative w/implicit electrodiffusion + ! 4 = Boussinesq w/implicit electrodiffusion + ! 5 = Inertial midpoint lowMach + ! 6 = Boussinesq inertial midpoint (electroneutral or electro-explicit) + + + ! Controls for deciding whether to plot various quantities (eg averaged, time averaged) + !---------------------- + plot_gradEpot = T + plot_umac_tavg = T + plot_Epot_tavg = T + plot_rho_tavg = F + + ! random number seed + ! 0 = unpredictable seed based on clock + ! positive = fixed seed + ! negative = (only for a restart) RNGs status is restored from checkpoint data + use_bl_rng = T + seed_momentum = 0 + seed_diffusion= 0 + seed_reaction = 0 + + ! Viscous friction L phi operator + ! if abs(visc_type) = 1, L = div beta grad + ! if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ] + ! if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ] + ! positive = assume constant coefficients + ! negative = assume spatially-varying coefficients + visc_type = 1 + visc_coef = 1 ! momentum diffusion coefficient 'eta' (dynamic viscosity), water=0.9d-2 + + advection_type = 0 ! 0 = centered explicit + ! 1 = unlimited bilinear bds in space and time + ! 2 = limited bliniear bds in space and time + ! 3 = unlimited quadratic bds in space and time + ! 4 = limited quadratic bds in space and time + + ! Stochastic momentum flux controls: + filtering_width = 0 ! If positive the *momentum* stochastic fluxes will be filtered (smoothed) + ! Stochastic *mass* fluxes are not filtered + stoch_stress_form = 1 ! 0=nonsymmetric (div(v)=0), 1=symmetric (no bulk) + + ! Initial conditions + !---------------------- + u_init(1:2) = 0.d0 0.d0 ! controls initial velocity + smoothing_width = 0.d0 ! scale factor for smoothing initial profile + initial_variance_mom = 0.d0 ! multiplicative factor for initial fluctuations + ! (if negative, total momentum is set to zero) + + ! Boundary conditions + !---------------------- + ! BC specifications: + ! -1 = periodic + ! 100 = no-slip wall (Dir condition for normal vel; Dir velocity condition for trans vel) + ! 101 = no-slip reservoir (Dir condition for normal vel; Dir velocity condition for trans vel) + ! 200 = slip wall (Dir condition for normal vel; Dir traction condition for trans vel) + ! 201 = slip reservoir (Dir condition for normal vel; Dir traction condition for trans vel) + ! For a complete list see bc.f90 + bc_lo(1:2) = -1 100 ! periodic on left, no slip on bottom + bc_hi(1:2) = -1 100 ! periodic on right, no slip on top + + ! Control for analyze_spectra.90 for calling HydroGrid + !---------------------- + hydro_grid_int = 0 ! How often to call updateHydroGrid + ! 0 if never + ! negative for projectHydroGrid custom analysis + ! positive for updateHydroGrid + + project_dir = 2 ! Projection direction (1=x, 2=y, 3=z) + ! Meaning: 0=analyze 3D data only (no projection needed for HydroGrid, + ! but still need projection if stats_int>0) + ! +dim=project along dim then analyze 2D only, + ! -dim=analyze 3D and then project along dim so we also analyze 2D data + ! It is better to use the conserved variables but it does not quite work for staggered + + !max_grid_projection(1:1) = 1024 ! parallelization box size for grid projection + max_grid_projection(1:1) = 32 ! parallelization box size for grid projection + stats_int = 10 ! Project grid for analysis + ! If positive, how often to compute mean and + ! standard deviation over reduced dimensions + stat_save_type = 0 ! Save instantaneous fields if = 0, otherwise time averages + + n_steps_save_stats = -1 ! How often to dump hstat/vstat and/or HydroGrid output files + n_steps_skip = 0 ! How many steps to skip at the beginning + analyze_conserved = F ! Should we use conserved variables for the analysis + ! (does not work well) + center_snapshots = F ! Should we use cell-centered momenta for the analysis + ! (will smooth fluctuations) + +/ + +&probin_chemistry + +/ + +&probin_multispecies + + ! Physical properties: + !---------------------- + fraction_tolerance = 1.d-12 ! For roundoff errors in mass and mole fractions + start_time = 0.d0 + + inverse_type = 1 ! Only for LAPACK: 1=inverse, 2=pseudo inverse + correct_flux = T ! Manually ensure mass is conserved to roundoff + print_error_norms = T + is_ideal_mixture = T ! If T assume Gamma=I (H=0) and simplify + is_nonisothermal = F ! If T Soret effect will be included + use_lapack = F ! Use LAPACK or iterative method for diffusion matrix (recommend False) + chi_iterations = 10 ! number of iterations used in Dbar2chi_iterative + + ! Initial and boundary conditions + !---------------------- + T_init(1:2) = 1.0 1.0 ! initial values for temperature (bottom/top, inside/outside circle, etc.) + temp_type = 0 ! for initializing temperature + + ! initial values for CONCENTRATIONS + c_init(1,1:3) = 0.005 0.005 0.99 + c_init(2,1:3) = 0.005 0.005 0.99 + + ! Thermodynamic and transport properties: + !---------------------- + + ! These are lower-triangules of symmetric matrices represented as vectors + ! Number of elements is (nspecies*(nspecies-1)/2) + ! The values are read row by row starting from top going down (this allows easy addition/deletion of new species/rows) + ! So D_12; D_13, D_23; D_14, D_24, D_34; ... + ! At low dilution (classical PNP), with last species being solvent + ! The last row are the self diffusion coefficients of the dilute solutes + ! Other entries are set to: D_ij = D_i*D_j / D_solvent + Dbar(1:3) = 1.33333333, 1.0 3.0 ! Just an example: D_water=1.5, D_ion1=1, D_ion2=2, D_12=1*2/1.5=4/3 + + ! Other examples: + + ! species = (1=Na, 2=Cl, 3=H2O) + ! So D_12; D_13, D_23; D_14, D_24, D_34; ... + ! diffusion of Na: 13.3e-6 ; Cl: 20.3e-6; self-water: 2.3e-5 + !Dbar(1:3) = D_12=1.17d-5, 1.33d-5 2.03d-5 + + ! species = (1=HCl, 2=NaOH, 3=NaCl, 4=H2O) + ! Binary ambipolar self-diffusion: HCl=3.336000e-05 NaOH=2.129000e-05 NaCl=1.611000e-05 + ! self-diffusion water = 2.3e-5 + !Dbar(1:6) = D_12=3.087976e-05, D_13=2.336650e-05 D_23=1.491226e-05, D_14=3.336000e-05 D_24=2.129000e-05 D_34=1.611000e-05 + + ! species = (1=Na+, 2=Cl-, 3=H+, 4=OH-, 5=H2O) + ! Self-diffusion: Na+ = 1.33e-5, Cl = 2.03e-5, H+ = 9.35e-5, OH = 5.33e-5, H20=2.3e-5 + !Dbar(1:10) = 0.1173869565e-4, 0.5408879652e-4 0.8255658417e-4, 0.3083589699e-4 0.4706531648e-4 0.2168644968e-3, 0.133e-4, 0.203e-4, 0.9353701657e-4, 0.5332523540e-4 + + plot_stag = F ! plot staggered velocities in separate plotfile + +/ + +&probin_charged + + use_charged_fluid = T + print_debye_len = T + dielectric_const = 1 ! epsilon (absolute units) water=708.01d-21 + charge_per_mass(1:3) = 1.0d0 -1.0d0 0.0d0 ! charge of each species + + ! 1 = Dirichlet (fixed potential) + ! 2 = Neumann (fixed charge density) + Epot_wall_bc_type(1,2) = 2 + Epot_wall_bc_type(2,2) = 2 + Epot_wall(1:2,2) = 1.d0 -1.d0 ! Potential or charge density at boundary + + ! Applied field along x: + E_ext_type = 1 + E_ext_value(1:2) = -1.0d0 0.d0 + + ! Notes from Andy Nonaka for Neumann: + ! The sign convention is a bit wonky as it represents dphi/dn (n is inward normal so that's why they are always the same sign). + ! The magnitude you want for solvability is equal to q^tot / (epsilon * L^tot) where L^tot is the total surface area of top and bottom walls. + ! q^tot can be found from the report at the beginning of standard out - doing this calculation analytically gives you results that are off in the 10th digit or so which is often enough to prevent the system from being solvable. + + electroneutral = F ! Use electroneutral formulation or PNP? + + epot_mg_verbose = 0 + epot_mg_abs_tol = 0 + epot_mg_rel_tol = 1.d-6 + + +/ + +! Stokes solver for velocity +&probin_gmres + + ! preconditioner type + ! 1 = projection preconditioner + !-1 = projection preconditioner with expensive pressure update + ! 2 = lower triangular preconditioner + !-2 = lower triangular preconditioner with negative sign + ! 3 = upper triangular preconditioner + !-3 = upper triangular preconditioner with negative sign + ! 4 = Block diagonal preconditioner + !-4 = Block diagonal preconditioner with negative sign + precon_type = 1 + + ! weighting of pressure when computing norms and inner products + p_norm_weight = 1.d0 + + ! scale theta_alpha, beta, gamma, and b_u by this, and then scale x_p by the inverse + scale_factor = 1.d0 + + ! MAC projection solver parameters: + mg_verbose = 0 ! multigrid verbosity + cg_verbose = 0 ! BiCGStab (mg_bottom_solver=1) verbosity + mg_max_vcycles = 1 ! maximum number of V-cycles + mg_minwidth = 2 ! length of box at coarsest multigrid level + mg_bottom_solver = 4 ! bottom solver type + ! 0 = smooths only, controlled by mg_nsmooths_bottom + ! 1 = BiCGStab + ! 4 = Fancy bottom solve that coarsens down additionally + ! and then applies mg_nsmooths_bottom smooths + mg_nsmooths_down = 2 ! number of smooths at each level on the way down + mg_nsmooths_up = 2 ! number of smooths at each level on the way up + mg_nsmooths_bottom = 8 ! number of smooths at the bottom (only if mg_bottom_solver=0) + mg_max_bottom_nlevels = 10 ! for mg_bottom_solver=4, number of additional levels of multigrid + mg_rel_tol = 1.d-6 ! relative tolerance stopping criteria + + ! Staggered multigrid solver parameters + stag_mg_verbosity = 0 ! verbosity (set to 1 for debugging the gmres solver!) + stag_mg_max_vcycles = 1 ! max number of v-cycles + stag_mg_minwidth = 2 ! length of box at coarsest multigrid level + stag_mg_bottom_solver = 4 ! bottom solver type + ! 0 = smooths only, controlled by mg_nsmooths_bottom + ! 4 = Fancy bottom solve that coarsens additionally + ! and then applies stag_mg_nsmooths_bottom smooths + stag_mg_nsmooths_down = 2 ! number of smooths at each level on the way down + stag_mg_nsmooths_up = 2 ! number of smooths at each level on the way up + stag_mg_nsmooths_bottom = 8 ! number of smooths at the bottom + stag_mg_max_bottom_nlevels = 10 ! for stag_mg_bottom_solver=4, number of additional levels of multigrid + stag_mg_omega = 1.d0 ! weighted-jacobi omega coefficient + stag_mg_smoother = 1 ! 0 = jacobi; 1 = 2*dm-color Gauss-Seidel + stag_mg_rel_tol = 1.d-6 ! relative tolerance stopping criteria + + ! GMRES solver parameters + gmres_rel_tol = 1.d-4 ! relative tolerance stopping criteria + gmres_abs_tol = 0.d0 ! absolute tolerance stopping criteria + gmres_verbose = 0 ! gmres verbosity; if greater than 1, more residuals will be printed out + gmres_max_outer = 20 ! max number of outer iterations + gmres_max_inner = 5 ! max number of inner iterations, or restart number + gmres_max_iter = 50 ! max number of gmres iterations + gmres_min_iter = 1 ! min number of gmres iterations + +/ diff --git a/staggered_grid/exec/chargedFluid/test/inputs_galen b/staggered_grid/exec/chargedFluid/test/inputs_galen new file mode 100644 index 00000000..c18c7b2f --- /dev/null +++ b/staggered_grid/exec/chargedFluid/test/inputs_galen @@ -0,0 +1,287 @@ +! debye ~ 8.2e-8 +! time step contraints +! electro: 3.5e-10 +! mass diff: 5.6e-12 +! advective: 2e-9? (based on max_vel=10) + + +! WARNING: gfortran has shown some issues reading matrices of rank higher than 2 from namelists +! be careful and check values are read correctly (e.g., c_bc) + +! Options for low Mach equations +&probin_common + + ! Problem specification + !---------------------- + dim_in = 2 ! 2D or 3D + prob_lo(1:3) = 0.d0 0.d0 0.d0 ! physical lo coordinate + prob_hi(1:3) = 2.636d-6 6.59d-7 2.059375d-8 ! physical hi coordinate + + n_cells(1:2) = 128 32 ! number of cells in domain + max_grid_size(1:2) = 32 32 ! max number of cells in a box + fixed_dt = 5.d-12 ! time step (if positive, fixed) + max_step = 10000 ! maximum number of time steps + plot_int = 100 ! Interval for writing a plotfile (for visit/amrvis) + + n_cells(1:2) = 256 64 ! number of cells in domain + max_grid_size(1:2) = 64 64 ! max number of cells in a box + fixed_dt = 1.25d-12 ! time step (if positive, fixed) + max_step = 50000 ! maximum number of time steps + plot_int = 500 ! Interval for writing a plotfile (for visit/amrvis) + + n_cells(1:2) = 512 128 ! number of cells in domain + max_grid_size(1:2) = 128 128 ! max number of cells in a box + fixed_dt = 3.125d-13 ! time step (if positive, fixed) + max_step = 250000 ! maximum number of time steps + plot_int = 2500 ! Interval for writing a plotfile (for visit/amrvis) + + ! Time-step control + !---------------------- + + ! Controls for number of steps between actions (for HydroGrid see below) + !---------------------- + chk_int = -1 ! Interval for writing a checkpoint + prob_type = 2 ! sets scalars, m, coefficients (see init.f90) + restart = -1 ! checkpoint restart number + print_int = 10 ! how often to output EOS drift and sum of conserved quantities + project_eos_int = 1000 ! how often to call project_onto_eos + + ! Physical parameters + !-------------------- + grav(1:2) = 0.d0 0.d0 ! gravity vector (negative is downwards) + molmass(1:3) = 5.887d-23 5.887d-23 2.992d-23 + rhobar(1:3) = 1.d0 1.d0 1.d0 + + ! stochastic forcing amplitudes (1 for physical values, 0 to run them off) + variance_coef_mom = 0.d0 + variance_coef_mass = 0.d0 + k_B = 1.38064852d-16 ! Boltzmann's constant + + + ! Algorithm control / selection + !---------------------- + algorithm_type = 0 ! differs from code to code. In binary and multispcies code: + ! 0 = Inertial algorithm + ! 2 = Overdamped with 2 RNGs + ! 3 = Iterative implicit + + ! Controls for deciding whether to plot various quantities (eg averaged, time averaged) + !---------------------- + plot_gradEpot = T + plot_umac_tavg = T + plot_Epot_tavg = T + plot_rho_tavg = T + + ! random number seed + ! 0 = unpredictable seed based on clock + ! positive = fixed seed + ! negative = (only for a restart) RNGs status is restored from checkpoint data + use_bl_rng = T + seed_momentum = 1 + seed_diffusion= 1 + seed_reaction = 1 + + ! Viscous friction L phi operator + ! if abs(visc_type) = 1, L = div beta grad + ! if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ] + ! if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ] + ! positive = assume constant coefficients + ! negative = assume spatially-varying coefficients + visc_type = 1 + visc_coef = 0.9d-2 ! momentum diffusion coefficient 'eta' + + advection_type = 0 ! 0 = centered explicit + ! 1 = unlimited bilinear bds in space and time + ! 2 = limited bliniear bds in space and time + ! 3 = unlimited quadratic bds in space and time + ! 4 = limited quadratic bds in space and time + + ! Stochastic momentum flux controls: + filtering_width = 0 ! If positive the *momentum* stochastic fluxes will be filtered (smoothed) + ! Stochastic *mass* fluxes are not filtered + stoch_stress_form = 1 ! 0=nonsymmetric (div(v)=0), 1=symmetric (no bulk) + + ! Initial conditions + !---------------------- + u_init(1:2) = 0.d0 0.d0 ! controls initial velocity + smoothing_width = 0.d0 ! scale factor for smoothing initial profile + initial_variance_mom = 0.d0 ! multiplicative factor for initial fluctuations + ! (if negative, total momentum is set to zero) + + ! Boundary conditions + !---------------------- + ! BC specifications: + ! -1 = periodic + ! 100 = no-slip wall (Dir condition for normal vel; Dir velocity condition for trans vel) + ! 101 = no-slip reservoir (Dir condition for normal vel; Dir velocity condition for trans vel) + ! 200 = slip wall (Dir condition for normal vel; Dir traction condition for trans vel) + ! 201 = slip reservoir (Dir condition for normal vel; Dir traction condition for trans vel) + ! For a complete list see bc.f90 + bc_lo(1:2) = -1 100 ! periodic on left, no slip on bottom + bc_hi(1:2) = -1 100 ! periodic on right, no slip on top + + ! Control for analyze_spectra.90 for calling HydroGrid + !---------------------- + hydro_grid_int = 0 ! How often to call updateHydroGrid + ! 0 if never + ! negative for projectHydroGrid custom analysis + ! positive for updateHydroGrid + + project_dir = 0 ! Projection direction (1=x, 2=y, 3=z) + ! Meaning: 0=analyze 3D data only (no projection needed for HydroGrid, + ! but still need projection if stats_int>0) + ! +dim=project along dim then analyze 2D only, + ! -dim=analyze 3D and then project along dim so we also analyze 2D data + ! It is better to use the conserved variables but it does not quite work for staggered + + max_grid_projection(1:1) = 128 ! parallelization parameters + stats_int = -1 ! Project grid for analysis + ! If positive, how often to compute mean and + ! standard deviation over reduced dimensions + n_steps_save_stats = -1 ! How often to dump HydroGrid output files + n_steps_skip = 0 ! How many steps to skip + analyze_conserved = F ! Should we use conserved variables for the analysis + ! (does not work well) + center_snapshots = F ! Should we use cell-centered momenta for the analysis + ! (will smooth fluctuations) + + nspecies = 3 +/ + +&probin_chemistry + +/ + +&probin_multispecies + + ! Physical properties: + !---------------------- + fraction_tolerance = 1.d-14 ! For roundoff errors in mass and mole fractions + start_time = 0.d0 + + inverse_type = 1 ! Only for LAPACK: 1=inverse, 2=pseudo inverse + correct_flux = T ! Manually ensure mass is conserved to roundoff + print_error_norms = T + is_ideal_mixture = T ! If T assume Gamma=I (H=0) and simplify + is_nonisothermal = F ! If T Soret effect will be included + use_lapack = F ! Use LAPACK or iterative method for diffusion matrix (recommend False) + chi_iterations = 10 ! number of iterations used in Dbar2chi_iterative + + ! Initial and boundary conditions + !---------------------- + T_init(1:2) = 300.d0 300.d0 ! initial values for temperature (bottom/top, inside/outside circle, etc.) + temp_type = 0 ! for initializing temperature + + ! initial values for CONCENTRATIONS + c_init(1,1:3) = 0.005 0.005 0.99 + c_init(2,1:3) = 0.005 0.005 0.99 + + ! Thermodynamic and transport properties: + !---------------------- + + ! These are lower-triangules of symmetric matrices represented as vectors + ! Number of elements is (nspecies*(nspecies-1)/2) + ! The values are read row by row starting from top going down (this allows easy addition/deletion of new species/rows) + ! So D_12; D_13, D_23; D_14, D_24, D_34; ... + Dbar(1:3) = 1.55e-5 1.89e-5 1.89e-5 + + plot_stag = F ! plot staggered velocities in separate plotfile + +/ + +&probin_charged + + use_charged_fluid = T + print_debye_len = T + dielectric_const = 708.01d-21 + charge_per_mass(1:3) = 2720.d0 -2720.d0 0.d0 + + induced_charge_eo = T + ac_iceo = F + + ! 1 = Dirichlet (fixed potential) + ! 2 = Neumann (fixed charge density) + Epot_wall_bc_type(1,2) = 1 + Epot_wall_bc_type(2,2) = 2 + Epot_wall(1,2) = 0.d0 + Epot_wall(2,2) = 0.d0 + + ! zero out eps on portions of a Dirichlet face to mimic dphi/dn=0 + zero_eps_on_wall_type = 1 + zero_eps_on_wall_left_end = 0.40d0 ! eg if set to 0.25, then eps will be set to 0 on the wall from 0*Lx --> 0.25*Lx + zero_eps_on_wall_right_start = 0.60d0 ! eg if set to 0.75, then eps will be set to 0 on the wall from 0.75*Lx --> 1.*Lx + E_ext_type = 1 + E_ext_value(1:2) = -1.5d9 0.d0 + E_ext_value(1:2) = -1.5d10 0.d0 + E_ext_value(1:2) = -1.5d11 0.d0 + E_ext_value(1:2) = -7.5d11 0.d0 + E_ext_value(1:2) = -1.5d12 0.d0 + E_ext_value(1:2) = -2.5d12 0.d0 + E_ext_value(1:2) = -3.0d12 0.d0 + E_ext_value(1:2) = -5.0d12 0.d0 + E_ext_value(1:2) = -1.0d13 0.d0 + E_ext_value(1:2) = -1.5d13 0.d0 +/ + +! Stokes solver for velocity +&probin_gmres + + ! preconditioner type + ! 1 = projection preconditioner + !-1 = projection preconditioner with expensive pressure update + ! 2 = lower triangular preconditioner + !-2 = lower triangular preconditioner with negative sign + ! 3 = upper triangular preconditioner + !-3 = upper triangular preconditioner with negative sign + ! 4 = Block diagonal preconditioner + !-4 = Block diagonal preconditioner with negative sign + precon_type = 1 + + ! weighting of pressure when computing norms and inner products + p_norm_weight = 1.d0 + + ! scale theta_alpha, beta, gamma, and b_u by this, and then scale x_p by the inverse + scale_factor = 1.d0 + + ! MAC projection solver parameters: + mg_verbose = 0 ! multigrid verbosity + cg_verbose = 0 ! BiCGStab (mg_bottom_solver=1) verbosity + mg_max_vcycles = 1 ! maximum number of V-cycles + mg_minwidth = 2 ! length of box at coarsest multigrid level + mg_bottom_solver = 4 ! bottom solver type + ! 0 = smooths only, controlled by mg_nsmooths_bottom + ! 1 = BiCGStab + ! 4 = Fancy bottom solve that coarsens down additionally + ! and then applies mg_nsmooths_bottom smooths + mg_nsmooths_down = 2 ! number of smooths at each level on the way down + mg_nsmooths_up = 2 ! number of smooths at each level on the way up + mg_nsmooths_bottom = 8 ! number of smooths at the bottom (only if mg_bottom_solver=0) + mg_max_bottom_nlevels = 10 ! for mg_bottom_solver=4, number of additional levels of multigrid + mg_rel_tol = 1.d-10 ! relative tolerance stopping criteria + + ! Staggered multigrid solver parameters + stag_mg_verbosity = 0 ! verbosity (set to 1 for debugging the gmres solver!) + stag_mg_max_vcycles = 1 ! max number of v-cycles + stag_mg_minwidth = 2 ! length of box at coarsest multigrid level + stag_mg_bottom_solver = 4 ! bottom solver type + ! 0 = smooths only, controlled by mg_nsmooths_bottom + ! 4 = Fancy bottom solve that coarsens additionally + ! and then applies stag_mg_nsmooths_bottom smooths + stag_mg_nsmooths_down = 2 ! number of smooths at each level on the way down + stag_mg_nsmooths_up = 2 ! number of smooths at each level on the way up + stag_mg_nsmooths_bottom = 8 ! number of smooths at the bottom + stag_mg_max_bottom_nlevels = 10 ! for stag_mg_bottom_solver=4, number of additional levels of multigrid + stag_mg_omega = 1.d0 ! weighted-jacobi omega coefficient + stag_mg_smoother = 1 ! 0 = jacobi; 1 = 2*dm-color Gauss-Seidel + stag_mg_rel_tol = 1.d-10 ! relative tolerance stopping criteria + + ! GMRES solver parameters + gmres_rel_tol = 1.d-9 ! relative tolerance stopping criteria + gmres_abs_tol = 0.d0 ! absolute tolerance stopping criteria + gmres_verbose = 0 ! gmres verbosity; if greater than 1, more residuals will be printed out + gmres_max_outer = 20 ! max number of outer iterations + gmres_max_inner = 5 ! max number of inner iterations, or restart number + gmres_max_iter = 100 ! max number of gmres iterations + gmres_min_iter = 1 ! min number of gmres iterations + +/ diff --git a/staggered_grid/exec/reactDiff/test/GNUmakefile.Donev b/staggered_grid/exec/reactDiff/test/GNUmakefile.Donev index 0824a417..f7518b84 100644 --- a/staggered_grid/exec/reactDiff/test/GNUmakefile.Donev +++ b/staggered_grid/exec/reactDiff/test/GNUmakefile.Donev @@ -1,8 +1,9 @@ -# FBOXLIB_HOME defines the directory in which we will find all the amrex code -# If you set FBOXLIB_HOME as an environment variable, this line will be ignored -FBOXLIB_HOME := ${HPC_SRC}/FluctHydro/FBoxLib +# BOXLIB_HOME defines the directory in which we will find all the BoxLib code +# If you set BOXLIB_HOME as an environment variable, this line will be ignored +BOXLIB_HOME := ${HPC_SRC}/FluctHydro/BoxLib HYDROLIB_HOME := ${HPC_SRC}/FluctHydro/HydroGrid STAGGERED_HOME := ${HPC_SRC}/FluctHydro/FluctHydro/staggered_grid +REACTDIFF_HOME := $(STAGGERED_HOME)/exec/reactDiff/test NDEBUG := t MPI := t @@ -10,42 +11,35 @@ OMP := PROF := MKVERBOSE := t -#COMP := gfortran -#CCOMP := gcc -COMP := gfortran63 -CCOMP := gcc63 +COMP := gfortran61 +CCOMP := gcc61 # need this to compile bl_rng.f90 (module to store random number engines and distributions) CXX11 = t -include $(FBOXLIB_HOME)/Tools/F_mk/GMakedefs.mak +include $(BOXLIB_HOME)/Tools/F_mk/GMakedefs.mak -# For CIMS donev machine: -#xtr_libraries += -Wl,-rpath=/usr/local/pkg/gcc/4.8.2/lib64 # RedHat -xtr_libraries += -Wl,-rpath=/usr/local/stow/gcc-6.3.0/lib64 # CentOS +xtr_libraries += -Wl,-rpath=/usr/local/pkg/gcc/6.1.0/lib64 -include $(FBOXLIB_HOME)/Tools/F_mk/GMakedefs.mak +VPATH_LOCATIONS += $(REACTDIFF_HOME) include $(STAGGERED_HOME)/src_reactDiff/GPackage.mak VPATH_LOCATIONS += $(STAGGERED_HOME)/src_reactDiff -include $(STAGGERED_HOME)/src_chemistry/GPackage.mak -VPATH_LOCATIONS += $(STAGGERED_HOME)/src_chemistry - include $(STAGGERED_HOME)/src_common/GPackage.mak VPATH_LOCATIONS += $(STAGGERED_HOME)/src_common include $(HYDROLIB_HOME)/GPackage.mak VPATH_LOCATIONS += $(HYDROLIB_HOME) -include $(FBOXLIB_HOME)/Src/BaseLib/GPackage.mak -VPATH_LOCATIONS += $(FBOXLIB_HOME)/Src/BaseLib +include $(BOXLIB_HOME)/Src/LinearSolvers/F_MG/GPackage.mak +VPATH_LOCATIONS += $(BOXLIB_HOME)/Src/LinearSolvers/F_MG -include $(FBOXLIB_HOME)/Src/MultiGrid/GPackage.mak -VPATH_LOCATIONS += $(FBOXLIB_HOME)/Src/MultiGrid +include $(BOXLIB_HOME)/Src/F_BaseLib/GPackage.mak +VPATH_LOCATIONS += $(BOXLIB_HOME)/Src/F_BaseLib main.$(suf).exe: $(objects) $(LINK.f90) -o main.$(suf).exe $(objects) $(libraries) -include $(FBOXLIB_HOME)/Tools/F_mk/GMakerules.mak +include $(BOXLIB_HOME)/Tools/F_mk/GMakerules.mak diff --git a/staggered_grid/src_charged/electrodiffusive_mass_fluxdiv.f90 b/staggered_grid/src_charged/electrodiffusive_mass_fluxdiv.f90 index e0478f49..5e34641b 100644 --- a/staggered_grid/src_charged/electrodiffusive_mass_fluxdiv.f90 +++ b/staggered_grid/src_charged/electrodiffusive_mass_fluxdiv.f90 @@ -32,7 +32,7 @@ module electrodiffusive_mass_fluxdiv_module subroutine electrodiffusive_mass_fluxdiv(mla,rho,Temp,rhoWchi, & diff_mass_flux,diff_mass_fluxdiv, & stoch_mass_flux, & - dx,the_bc_tower,charge, & + dx,stage_time,the_bc_tower,charge, & grad_Epot,Epot,permittivity,dt, & zero_initial_Epot) @@ -47,6 +47,7 @@ subroutine electrodiffusive_mass_fluxdiv(mla,rho,Temp,rhoWchi, & type(multifab) , intent(inout) :: diff_mass_fluxdiv(:) type(multifab) , intent(in ) :: stoch_mass_flux(:,:) real(kind=dp_t), intent(in ) :: dx(:,:) + real(kind=dp_t), intent(in ) :: stage_time type(bc_tower) , intent(in ) :: the_bc_tower type(multifab) , intent(inout) :: charge(:) type(multifab) , intent(inout) :: grad_Epot(:,:) @@ -80,7 +81,7 @@ subroutine electrodiffusive_mass_fluxdiv(mla,rho,Temp,rhoWchi, & ! compute the face-centered electro_mass_flux (each direction: cells+1 faces while ! cells contain interior+2 ghost cells) call electrodiffusive_mass_flux(mla,rho,Temp,rhoWchi,electro_mass_flux, & - diff_mass_flux,stoch_mass_flux,dx,the_bc_tower, & + diff_mass_flux,stoch_mass_flux,dx,stage_time,the_bc_tower, & charge,grad_Epot,Epot,permittivity,dt,zero_initial_Epot) ! add fluxes to diff_mass_flux @@ -107,7 +108,7 @@ end subroutine electrodiffusive_mass_fluxdiv subroutine electrodiffusive_mass_flux(mla,rho,Temp,rhoWchi,electro_mass_flux, & diff_mass_flux,stoch_mass_flux, & - dx,the_bc_tower,charge,grad_Epot,Epot, & + dx,stage_time,the_bc_tower,charge,grad_Epot,Epot, & permittivity,dt,zero_initial_Epot) ! this computes "-F = A_Phi grad Phi" @@ -120,6 +121,7 @@ subroutine electrodiffusive_mass_flux(mla,rho,Temp,rhoWchi,electro_mass_flux, & type(multifab) , intent(in ) :: diff_mass_flux(:,:) type(multifab) , intent(in ) :: stoch_mass_flux(:,:) real(kind=dp_t), intent(in ) :: dx(:,:) + real(kind=dp_t), intent(in ) :: stage_time type(bc_tower) , intent(in ) :: the_bc_tower type(multifab) , intent(inout) :: charge(:) type(multifab) , intent(inout) :: grad_Epot(:,:) @@ -212,7 +214,7 @@ subroutine electrodiffusive_mass_flux(mla,rho,Temp,rhoWchi,electro_mass_flux, & ! fill ghost cells for Epot at walls using Dirichlet value call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in = stage_time, dx_in=dx(n,:)) ! set alpha=0 call setval(alpha(n),0.d0,all=.true.) @@ -361,7 +363,7 @@ subroutine electrodiffusive_mass_flux(mla,rho,Temp,rhoWchi,electro_mass_flux, & if (E_ext_type .ne. 0) then ! compute external electric field on edges - call compute_E_ext(mla,E_ext) + call compute_E_ext(mla,E_ext,stage_time,dt) ! compute epsilon*E_ext do n=1,nlevs @@ -448,7 +450,7 @@ subroutine electrodiffusive_mass_flux(mla,rho,Temp,rhoWchi,electro_mass_flux, & ! homogeneous BC's, this routine will properly fill the ghost cells do n=1,nlevs call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=stage_time, dx_in=dx(n,:)) call multifab_fill_boundary(Epot(n)) end do diff --git a/staggered_grid/src_charged/fluid_charge.f90 b/staggered_grid/src_charged/fluid_charge.f90 index fc066308..533743e4 100644 --- a/staggered_grid/src_charged/fluid_charge.f90 +++ b/staggered_grid/src_charged/fluid_charge.f90 @@ -6,6 +6,7 @@ module fluid_charge_module use div_and_grad_module use define_bc_module use bc_module + use inhomogeneous_bc_val_module use mass_flux_utilities_module use matvec_mul_module use compute_mixture_properties_module @@ -15,7 +16,7 @@ module fluid_charge_module use probin_charged_module, only: charge_per_mass, dpdt_factor, & dielectric_const, dielectric_type, & E_ext_type, E_ext_value, zero_eps_on_wall_type, & - zero_eps_on_wall_left_end, zero_eps_on_wall_right_start + zero_eps_on_wall_left_end, zero_eps_on_wall_right_start,ac_iceo implicit none @@ -959,10 +960,11 @@ end subroutine Lorentz_force_gradeps_3d end subroutine compute_Lorentz_force - subroutine compute_E_ext(mla,E_ext) + subroutine compute_E_ext(mla,E_ext,time,dt) type(ml_layout), intent(in ) :: mla type(multifab) , intent(inout) :: E_ext(:,:) + real(kind=dp_t), intent(in ) :: time,dt integer :: n, nlevs, dm, i @@ -973,7 +975,13 @@ subroutine compute_E_ext(mla,E_ext) do n=1,nlevs do i=1,dm - call multifab_setval(E_ext(n,i),E_ext_value(i),all=.true.) + ! if we're doing alternating current iceo, the applied external electric field + ! oscillates in time with amplitude = E_ext_value(i) + if (ac_iceo.and.(i.eq.1)) then + call multifab_setval(E_ext(n,i), alternating_current_efield(time),all=.true.) + else + call multifab_setval(E_ext(n,i),E_ext_value(i),all=.true.) + endif end do end do diff --git a/staggered_grid/src_charged/probin_charged.f90 b/staggered_grid/src_charged/probin_charged.f90 index 7fb292d9..7a109684 100644 --- a/staggered_grid/src_charged/probin_charged.f90 +++ b/staggered_grid/src_charged/probin_charged.f90 @@ -19,7 +19,7 @@ module probin_charged_module integer :: E_ext_type real(kind=dp_t) :: E_ext_value(1:3) logical :: electroneutral - logical :: induced_charge_eo + logical :: induced_charge_eo, ac_iceo integer :: zero_eps_on_wall_type integer :: zero_charge_on_wall_type real(kind=dp_t) :: zero_eps_on_wall_left_end, zero_eps_on_wall_right_start @@ -53,7 +53,8 @@ module probin_charged_module namelist /probin_charged/ E_ext_value ! spacedim-vector specifying external E field namelist /probin_charged/ electroneutral ! use electroneutral diffusion fluxes - namelist /probin_charged/ induced_charge_eo ! are we simulating ICEO? + namelist /probin_charged/ induced_charge_eo ! Are we simulating ICEO? + namelist /probin_charged/ ac_iceo ! If so, do we want to drive the flow with direct or alternating current? namelist /probin_charged/ relxn_param_charge ! Used to prevent slow buildup of charge for electroneutral, keep at 1.0 namelist /probin_charged/ zero_eps_on_wall_type ! set eps=0 on certain Dirichlet walls ! if we want homogeneous Neumann bc's on @@ -113,6 +114,7 @@ subroutine probin_charged_init() electroneutral = .false. induced_charge_eo = .false. + ac_iceo = .false. relxn_param_charge = 1.0d0 zero_eps_on_wall_type = 0 zero_charge_on_wall_type = 0 @@ -205,6 +207,11 @@ subroutine probin_charged_init() call get_command_argument(farg, value = fname) read(fname, *) induced_charge_eo + case ('--ac_iceo') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) ac_iceo + case ('--zero_eps_on_wall_type') farg = farg + 1 call get_command_argument(farg, value = fname) diff --git a/staggered_grid/src_common/analyze_spectra.f90 b/staggered_grid/src_common/analyze_spectra.f90 index bb9c954c..4010c26e 100644 --- a/staggered_grid/src_common/analyze_spectra.f90 +++ b/staggered_grid/src_common/analyze_spectra.f90 @@ -17,7 +17,7 @@ module analyze_spectra_module use utility_module use probin_common_module, only: n_cells, prob_lo, prob_hi, & hydro_grid_int, project_dir, max_grid_projection, stats_int, n_steps_save_stats, & - center_snapshots, analyze_conserved, histogram_unit, density_weights + center_snapshots, analyze_conserved, histogram_unit, density_weights,nspecies implicit none @@ -60,6 +60,7 @@ module analyze_spectra_module subroutine initialize_hydro_grid(mla,s_in,dt,dx,namelist_file, & nspecies_in, nscal_in, exclude_last_species_in, & analyze_velocity, analyze_density, analyze_temperature, & + analyze_charge_fluxes, & heat_capacity_in,structFactMultiplier) type(ml_layout), intent(in ) :: mla type(multifab) , intent(inout) :: s_in(:) ! A cell-centered multifab on the desired grid (to grab layout and grid size from) @@ -71,6 +72,7 @@ subroutine initialize_hydro_grid(mla,s_in,dt,dx,namelist_file, & integer , intent(in ) :: nscal_in ! Additional number of scalar fields if any logical, intent(in) :: exclude_last_species_in, analyze_velocity, analyze_density, analyze_temperature + logical, intent(in) :: analyze_charge_fluxes ! Pass exclude_last_species=.false. if you want to analyze all nspecies densities/concentrations ! Or pass exclude_last_species=.true. if you want to pass total density rho as the first scalar and ! then only include only n_species-1 additional partial densities/concentrations @@ -108,6 +110,9 @@ subroutine initialize_hydro_grid(mla,s_in,dt,dx,namelist_file, & if(analyze_temperature) nvar = nvar + 1 nscal_analysis = nvar ! Total number of scalar variables to analyze if(analyze_velocity) nvar = nvar + dm + if (analyze_charge_fluxes) then + nvar = nvar + dm*(nspecies+1) + endif if ( parallel_IOprocessor() ) then write(*,*) "Allocating storage for ", nvar, " variables to analyze using HydroGrid" end if @@ -381,17 +386,19 @@ subroutine analyze_hydro_grid(mla,dt,dx,step,umac,rho,temperature,scalars) ! This routines collects a bunch of different variables into a single multifab for easier analysis ! The components are ordered as, whenever present: ! umac, rho, rho_k (analyze_conserved=T) or c_k=rho_k/rho, T, scalars - subroutine gather_hydro_grid(mla,la_s,s_hydro,umac,rho,temperature,scalars, variable_names) + subroutine gather_hydro_grid(mla,la_s,s_hydro,umac,rho,temperature,scalars,variable_names,charge_fluxes) type(ml_layout), intent(in ) :: mla ! Layout of hydrodynamic arrays type(layout) , intent(in ) :: la_s ! Layout of s_hydro (changes depending on context) type(multifab) , intent(inout) :: s_hydro ! Output: collected grid data for analysis type(multifab) , intent(in), optional :: umac(:,:) type(multifab) , intent(in), dimension(:), optional :: rho, temperature, scalars + type(multifab) , intent(in), optional :: charge_fluxes(:,:) character(len=20), intent(out), optional :: variable_names(nvar) ! Local integer :: i, ii, jj, n, nlevs, dm, pdim, comp, n_passive_scals, species type(multifab) :: mac_cc(mla%nlevel) + type(multifab) :: charge_fluxes_cc(mla%nlevel) nlevs = mla%nlevel dm = mla%dim @@ -496,7 +503,52 @@ subroutine gather_hydro_grid(mla,la_s,s_hydro,umac,rho,temperature,scalars, vari comp=comp+1 n_passive_scals=n_passive_scals-1 end if - + + if (present(charge_fluxes)) then + ! build cell centered multifab for charge_Fluxes + do n=1,nlevs + call multifab_build(charge_fluxes_cc(n), mla%la(n), dm*(nspecies+1), 0) + end do + ! copy the data at the edges to the cell centers + do i=1,dm +! call shift_face_to_cc(mla,charge_fluxes(:,i),1,charge_fluxes_cc,i,nspecies+1) + ! AJN - I think you want this + call shift_face_to_cc(mla,charge_fluxes(:,i),1,charge_fluxes_cc,(i-1)*(nspecies+1)+1,nspecies+1) + enddo + ! copy cc data to s_hydro + call copy(s_hydro,comp,charge_fluxes_cc(1),1,dm*(nspecies+1)) + ! destroy the cc multifab + do n=1,nlevs + call multifab_destroy(charge_fluxes_cc(n)) + end do + ! write variable names + if (present(variable_names)) then + ! write variable names in x + do species=1,nspecies + write(variable_names(comp+species-1),"(A,I0)") "chrg_flx_x_", species + end do + write(variable_names(comp+nspecies),"(A,I0)") "tot_chrg_flx_x" + comp = comp + (nspecies+1) + + ! write variable names in y + do species=1,nspecies + write(variable_names(comp+species-1),"(A,I0)") "chrg_flx_y_", species + end do + write(variable_names(comp+nspecies),"(A,I0)") "tot_chrg_flx_y" + comp = comp + (nspecies+1) + + if (dm>2) then + ! write variable names in z + do species=1,nspecies + write(variable_names(comp+species-1),"(A,I0)") "chrg_flx_z_", species + end do + write(variable_names(comp+nspecies),"(A,I0)") "tot_chrg_flx_z" + comp = comp + (nspecies+1) + endif + endif + endif + + if(present(scalars)) then if(n_passive_scals<1) then call parallel_abort("No room left to store passive scalars, nscal_analysis too small") @@ -820,13 +872,14 @@ subroutine average_1D(variables_1D) ! ! Next, it writes out "horizontal" averages. ! Thus, in both 2D and 3D, it writes out a 1D text file called "hstatXXXXXX" - subroutine print_stats(mla,dx,step,time,umac,rho,temperature,scalars) + subroutine print_stats(mla,dx,step,time,umac,rho,temperature,scalars,charge_fluxes) type(ml_layout), intent(in ) :: mla real(dp_t) , intent(in ) :: dx(:,:) integer , intent(in ) :: step real(kind=dp_t), intent(in ) :: time type(multifab) , intent(in), optional :: umac(:,:) type(multifab) , intent(in), dimension(:), optional :: rho,temperature,scalars + type(multifab) , intent(in), optional :: charge_fluxes(:,:) ! Local integer nlevs,dm,pdim,qdim,qdim2,i,n,comp @@ -869,7 +922,8 @@ subroutine print_stats(mla,dx,step,time,umac,rho,temperature,scalars) ! Re-distribute the full grid into a grid of "tall skinny boxes" ! These boxes are not distributed along project_dim so we can do local analysis easily ! ------------------------- - call gather_hydro_grid(mla,la_dir,s_dir,umac,rho,temperature,scalars, variable_names(1:nvar)) + call gather_hydro_grid(mla,la_dir,s_dir,umac,rho,temperature,scalars, & + variable_names(1:nvar),charge_fluxes) ! Compute s_projected (average) and s_var (variance) ! ------------------------- diff --git a/staggered_grid/src_common/probin_common.f90 b/staggered_grid/src_common/probin_common.f90 index c8b99dd2..fa01880c 100644 --- a/staggered_grid/src_common/probin_common.f90 +++ b/staggered_grid/src_common/probin_common.f90 @@ -25,8 +25,9 @@ module probin_common_module real(dp_t),save :: wallspeed_lo(MAX_SPACEDIM-1,MAX_SPACEDIM) real(dp_t),save :: wallspeed_hi(MAX_SPACEDIM-1,MAX_SPACEDIM) integer,save :: hydro_grid_int,project_dir,max_grid_projection(2) - integer,save :: stats_int,n_steps_save_stats,n_steps_skip,histogram_unit,reset_tavg_step,stat_save_type - logical,save :: analyze_conserved,center_snapshots,reset_tavg_vals + integer,save :: stats_int,n_steps_save_stats,n_steps_skip,histogram_unit,stat_save_type + logical,save :: reset_averages + logical,save :: analyze_conserved,center_snapshots real(dp_t),save :: variance_coef_mom,variance_coef_mass,initial_variance_mom,initial_variance_mass real(dp_t),save :: k_B,Runiv,visc_coef integer,save :: stoch_stress_form,filtering_width,max_step @@ -38,6 +39,8 @@ module probin_common_module real(dp_t),save :: density_weights(0:MAX_SPECIES) integer,save :: shift_cc_to_boundary(MAX_SPACEDIM,2) logical,save :: plot_avg_gradPhiApprox, plot_gradEpot + logical,save :: plot_mass_fluxes, plot_mass_fluxes_tavg + logical,save :: plot_charge_fluxes, plot_charge_fluxes_tavg logical,save :: plot_averaged_vel, plot_shifted_vel logical,save :: plot_umac_tavg, plot_Epot_tavg, plot_rho_tavg logical,save :: plot_debug @@ -122,6 +125,10 @@ module probin_common_module namelist /probin_common/ plot_averaged_vel ! Write spatially-averaged velocities namelist /probin_common/ plot_avg_gradPhiApprox ! Write ambipolar approximation to gradPhi or not namelist /probin_common/ plot_gradEpot ! Write out gradient of electric potential + namelist /probin_common/ plot_mass_fluxes ! Write out mass fluxes of individual species, plus the total mass flux + namelist /probin_common/ plot_charge_fluxes ! Write out charge fluxes of individual species, plus the total charge flux + namelist /probin_common/ plot_mass_fluxes_tavg ! Write out mass fluxes of individual species, plus the total mass flux, time_averaged + namelist /probin_common/ plot_charge_fluxes_tavg ! Write out charge fluxes of individual species, plus the total charge flux, time_averaged namelist /probin_common/ plot_debug ! for boussinesq, print out rho. for electroneutral, print out charge @@ -220,8 +227,7 @@ module probin_common_module ! it's set to 0, we save instantaneous fields; else, we save time averaged fields namelist /probin_common/ n_steps_save_stats ! How often to dump HydroGrid output files namelist /probin_common/ n_steps_skip ! How many steps to skip - namelist /probin_common/ reset_tavg_vals ! Boolean to indicate if we want to reset the time averaged vals or not - namelist /probin_common/ reset_tavg_step ! If reset_tavg_vals = T, at what timestep do we reset + namelist /probin_common/ reset_averages ! reset time-averaged quantities every n_steps_skip namelist /probin_common/ analyze_conserved ! Should we use conserved variables for the analysis ! (does not work well) namelist /probin_common/ center_snapshots ! Should we use cell-centered momenta for the analysis @@ -310,6 +316,10 @@ subroutine probin_common_init() plot_avg_gradPhiApprox = .false. plot_gradEpot = .false. plot_debug = .false. + plot_mass_fluxes = .false. + plot_charge_fluxes = .false. + plot_mass_fluxes_tavg = .false. + plot_charge_fluxes_tavg = .false. barodiffusion_type = 0 @@ -349,8 +359,7 @@ subroutine probin_common_init() stat_save_type = 0 ! by default, vstat/hstat files contain instantaneous fields, not time-averaged ones. n_steps_save_stats = -1 n_steps_skip = 0 - reset_tavg_vals = .false. - reset_tavg_step = 0 + reset_averages = .false. analyze_conserved = .false. center_snapshots = .false. analyze_cuts = 0 @@ -614,6 +623,26 @@ subroutine probin_common_init() call get_command_argument(farg, value = fname) read(fname, *) plot_debug + case ('--plot_mass_fluxes') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) plot_mass_fluxes + + case ('--plot_charge_fluxes') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) plot_charge_fluxes + + case ('--plot_mass_fluxes_tavg') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) plot_mass_fluxes_tavg + + case ('--plot_charge_fluxes_tavg') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) plot_charge_fluxes_tavg + case ('--barodiffusion_type') farg = farg + 1 call get_command_argument(farg, value = fname) @@ -823,6 +852,11 @@ subroutine probin_common_init() call get_command_argument(farg, value = fname) read(fname, *) n_steps_skip + case ('--reset_averages') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) reset_averages + case ('--analyze_conserved') farg = farg + 1 call get_command_argument(farg, value = fname) diff --git a/staggered_grid/src_energy/main_driver.f90 b/staggered_grid/src_energy/main_driver.f90 index 6ab99932..acb778f6 100644 --- a/staggered_grid/src_energy/main_driver.f90 +++ b/staggered_grid/src_energy/main_driver.f90 @@ -382,7 +382,8 @@ subroutine main_driver() exclude_last_species_in=.false., & analyze_velocity=.true., & analyze_density=.true., & - analyze_temperature=.true.) + analyze_temperature=.true., & + analyze_Epot=.false.) close(unit=un) end if diff --git a/staggered_grid/src_lowMach/advance_timestep_imp_bousq.f90 b/staggered_grid/src_lowMach/advance_timestep_imp_bousq.f90 index 4e67e57d..cf7fd32e 100644 --- a/staggered_grid/src_lowMach/advance_timestep_imp_bousq.f90 +++ b/staggered_grid/src_lowMach/advance_timestep_imp_bousq.f90 @@ -307,7 +307,7 @@ subroutine advance_timestep_imp_bousq(mla,umac,rho_old,rho_new,rhotot_old,rhotot call multifab_setval(Epot(n),0.d0,all=.true.) ! fill ghost cells for Epot at walls using Dirichlet value call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=time, dx_in=dx(n,:)) end do ! for inhomogeneous Neumann bc's for electric potential, put in homogeneous form @@ -344,7 +344,7 @@ subroutine advance_timestep_imp_bousq(mla,umac,rho_old,rho_new,rhotot_old,rhotot call multifab_sub_sub_s(Epot(1),sum) do n=1,nlevs call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=time, dx_in=dx(n,:)) call multifab_fill_boundary(Epot(n)) end do end if @@ -663,7 +663,7 @@ subroutine advance_timestep_imp_bousq(mla,umac,rho_old,rho_new,rhotot_old,rhotot call multifab_setval(Epot(n),0.d0,all=.true.) ! fill ghost cells for Epot at walls using Dirichlet value call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=time, dx_in=dx(n,:)) end do ! for inhomogeneous Neumann bc's for electric potential, put in homogeneous form @@ -700,7 +700,7 @@ subroutine advance_timestep_imp_bousq(mla,umac,rho_old,rho_new,rhotot_old,rhotot call multifab_sub_sub_s(Epot(1),sum) do n=1,nlevs call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=time, dx_in=dx(n,:)) call multifab_fill_boundary(Epot(n)) end do end if diff --git a/staggered_grid/src_lowMach/advance_timestep_inertial.f90 b/staggered_grid/src_lowMach/advance_timestep_inertial.f90 index 5b551cf9..446eec66 100644 --- a/staggered_grid/src_lowMach/advance_timestep_inertial.f90 +++ b/staggered_grid/src_lowMach/advance_timestep_inertial.f90 @@ -45,6 +45,7 @@ module advance_timestep_inertial_module subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_new, & gradp_baro,pi,eta,eta_ed,kappa,Temp,Temp_ed, & diff_mass_fluxdiv,stoch_mass_fluxdiv,stoch_mass_flux, & + total_mass_flux_w_adv, & dx,dt,time,the_bc_tower,istep, & grad_Epot_old,grad_Epot_new,charge_old,charge_new, & Epot,permittivity) @@ -66,6 +67,8 @@ subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_ type(multifab) , intent(inout) :: diff_mass_fluxdiv(:) type(multifab) , intent(inout) :: stoch_mass_fluxdiv(:) type(multifab) , intent(inout) :: stoch_mass_flux(:,:) + type(multifab) , intent(inout) :: total_mass_flux_w_adv(:,:) ! the actual total mass fluxes, including advective part. different than the + ! (somewhat confusingly named) total_mass_flux defined locally below real(kind=dp_t), intent(in ) :: dx(:,:),dt,time type(bc_tower) , intent(in ) :: the_bc_tower integer , intent(in ) :: istep @@ -103,11 +106,14 @@ subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_ ! only used when variance_coef_mom>0 type(multifab) :: stoch_mom_fluxdiv(mla%nlevel,mla%dim) + ! only used when advection_type = 0 ie centered explicit + type(multifab) :: adv_mass_flux(mla%nlevel,mla%dim) !SC + ! only used when use_charged_fluid=T type(multifab) :: Lorentz_force_old(mla%nlevel,mla%dim) type(multifab) :: Lorentz_force_new(mla%nlevel,mla%dim) - integer :: i,dm,n,nlevs + integer :: i,dm,n,nlevs,comp real(kind=dp_t) :: theta_alpha, norm_pre_rhs, gmres_abs_tol_in @@ -154,6 +160,14 @@ subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_ end do end do end if + + if (advection_type.eq.0) then !SC + do n=1,nlevs + do i=1,dm + call multifab_build_edge( adv_mass_flux(n,i),mla%la(n),nspecies,0,i) + enddo + enddo + endif if (use_charged_fluid) then do n=1,nlevs @@ -374,7 +388,7 @@ subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_ end if end do end do - + ! set the Dirichlet velocity value on reservoir faces call reservoir_bc_fill(mla,total_mass_flux,vel_bc_n,the_bc_tower%bc_tower_array) @@ -734,6 +748,54 @@ subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_ end do end do + + ! compute advective mass flux to be added to total_mass_flux + if (advection_type.eq.0) then + do n=1,nlevs + do i=1,dm + do comp=1,nspecies + call multifab_copy_c(adv_mass_flux(n,i),comp,umac(n,i),1,1,0) + end do + call multifab_mult_mult_c(adv_mass_flux(n,i),1,rho_fc(n,i),1,nspecies,0) + enddo + enddo + + ! Now say (total_mass_flux_w_adv = total_mass_flux - advective_mass_flux) + ! + ! total_mass_flux_w_adv = +[stuff]*grad(rho) + [stuff]*grad(epot) + [stochastic flux] - rho_tot*v + ! + ! BUT we flip the sign below, see the comment there. + do n=1,nlevs + do i=1,dm + ! first copy [diffusive_flux] + [stoch_flux] from above to total_mass_flux_w_adv + call multifab_copy_c(total_mass_flux_w_adv(n,i),1,total_mass_flux(n,i),1,nspecies,0) + + ! now subtract off the advective flux from the diffusive and stochastic fluxes to get total mass flux + call multifab_sub_sub_c(total_mass_flux_w_adv(n,i),1,adv_mass_flux(n,i),1,nspecies,0) + end do + end do + + ! now that we have all the individual species mass fluxes--add them all up to get the total mass flux + do n=1,nlevs + do i=1,dm + call multifab_copy_c(total_mass_flux_w_adv(n,i),nspecies+1,total_mass_flux_w_adv(n,i),1,1,0) + do comp=2,nspecies + call multifab_plus_plus_c(total_mass_flux_w_adv(n,i),nspecies+1,total_mass_flux_w_adv(n,i),comp,1,0) + enddo + enddo + enddo + + ! multiply the total mass flux w adv by -1, so that we take the convention that + ! the mass flux has the same sign as the fluid flow + do n=1,nlevs + do i=1,dm + call multifab_mult_mult_s_c(total_mass_flux_w_adv(n,i),1,-1.d0,nspecies+1,0) + enddo + enddo + + endif + + ! set the Dirichlet velocity value on reservoir faces call reservoir_bc_fill(mla,total_mass_flux,vel_bc_n,the_bc_tower%bc_tower_array) @@ -909,6 +971,14 @@ subroutine advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_ end do end if + if (advection_type.eq.0) then + do n=1,nlevs + do i=1,dm + call multifab_destroy(adv_mass_flux(n,i)) + end do + end do + end if + if (use_charged_fluid) then do n=1,nlevs do i=1,dm diff --git a/staggered_grid/src_lowMach/advance_timestep_iterative.f90 b/staggered_grid/src_lowMach/advance_timestep_iterative.f90 index 98f58793..7d2afc07 100644 --- a/staggered_grid/src_lowMach/advance_timestep_iterative.f90 +++ b/staggered_grid/src_lowMach/advance_timestep_iterative.f90 @@ -390,7 +390,7 @@ subroutine advance_timestep_iterative(mla,umac,rho_old,rho_new,rhotot_old,rhotot call multifab_setval(Epot(n),0.d0,all=.true.) ! fill ghost cells for Epot at walls using Dirichlet value call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=time, dx_in=dx(n,:)) end do ! for inhomogeneous Neumann bc's for electric potential, put in homogeneous form @@ -427,7 +427,7 @@ subroutine advance_timestep_iterative(mla,umac,rho_old,rho_new,rhotot_old,rhotot call multifab_sub_sub_s(Epot(1),sum) do n=1,nlevs call multifab_physbc(Epot(n),1,Epot_bc_comp,1,the_bc_tower%bc_tower_array(n), & - dx_in=dx(n,:)) + time_in=time, dx_in=dx(n,:)) call multifab_fill_boundary(Epot(n)) end do end if diff --git a/staggered_grid/src_lowMach/checkpoint.f90 b/staggered_grid/src_lowMach/checkpoint.f90 index 403dfed8..ea03eed3 100644 --- a/staggered_grid/src_lowMach/checkpoint.f90 +++ b/staggered_grid/src_lowMach/checkpoint.f90 @@ -11,7 +11,7 @@ module checkpoint_module use fabio_module use probin_common_module, only: dim_in, use_bl_rng, nspecies, & seed_diffusion, seed_momentum, seed_reaction, & - check_base_name + check_base_name, advection_type use probin_chemistry_module, only: nreactions use probin_charged_module, only: use_charged_fluid @@ -22,7 +22,8 @@ module checkpoint_module public :: checkpoint_write, checkpoint_read contains - subroutine checkpoint_write(mla,rho,rho_sum,rhotot,pi,umac,umac_sum,Epot,Epot_sum,grad_Epot,time,dt,istep_to_write) + subroutine checkpoint_write(mla,rho,rho_sum,rhotot,pi,umac,umac_sum,mass_fluxes,mass_fluxes_sum, & + Epot,Epot_sum,grad_Epot,time,dt,istep_to_write) type(ml_layout), intent(in ) :: mla type(multifab) , intent(in ) :: rho(:) ! cell-centered partial densities @@ -34,6 +35,8 @@ subroutine checkpoint_write(mla,rho,rho_sum,rhotot,pi,umac,umac_sum,Epot,Epot_su type(multifab) , intent(in ) :: grad_Epot(:,:) ! edge-based electric force type(multifab) , intent(in ) :: umac(:,:) ! edge-based velocities type(multifab) , intent(in ) :: umac_sum(:,:) ! edge-based sum of velocities + type(multifab) , intent(in ) :: mass_fluxes(:,:) ! edge-based mass fluxes of each species + type(multifab) , intent(in ) :: mass_fluxes_sum(:,:) ! edge-based sum of mass fluxes of each species integer , intent(in ) :: istep_to_write real(kind=dp_t), intent(in ) :: time,dt @@ -116,21 +119,48 @@ subroutine checkpoint_write(mla,rho,rho_sum,rhotot,pi,umac,umac_sum,Epot,Epot_su ! staggered quantities (normal velocity) + ! only write the mass fluxes if advection type = 0 do n=1,nlevs do i=1,dm if (use_charged_fluid) then - call multifab_build_edge(chkdata_edge(n,i),mla%la(n),3,0,i) - call multifab_copy_c(chkdata_edge(n,i),1,umac(n,i),1,1) - call multifab_copy_c(chkdata_edge(n,i),2,umac_sum(n,i),1,1) - call multifab_copy_c(chkdata_edge(n,i),3,grad_Epot(n,i),1,1) + if (advection_type.eq.0) then + ! + ! umac 1 + ! umac_sum 1 + ! mass_fluxes nspecies+1 <--mass fluxes contain a comp for each species flux, plus one comp for their sum + ! mass_fluxes_sum nspecies+1 + ! grad_Epot 1 + call multifab_build_edge(chkdata_edge(n,i),mla%la(n),(nspecies+1)*2+3,0,i) + call multifab_copy_c(chkdata_edge(n,i),1,umac(n,i),1,1) + call multifab_copy_c(chkdata_edge(n,i),2,umac_sum(n,i),1,1) + !call multifab_copy_c(chkdata_edge(n,i),3,mass_fluxes(n,i),nspecies+1,1) + call multifab_copy_c(chkdata_edge(n,i),3,mass_fluxes(n,i),1,nspecies+1) + call multifab_copy_c(chkdata_edge(n,i),(nspecies+1)+3,mass_fluxes_sum(n,i),1,nspecies+1) + call multifab_copy_c(chkdata_edge(n,i),2*(nspecies+1)+3,grad_Epot(n,i),1,1) + + else + call multifab_build_edge(chkdata_edge(n,i),mla%la(n),3,0,i) + call multifab_copy_c(chkdata_edge(n,i),1,umac(n,i),1,1) + call multifab_copy_c(chkdata_edge(n,i),2,umac_sum(n,i),1,1) + call multifab_copy_c(chkdata_edge(n,i),3,grad_Epot(n,i),1,1) + endif else - call multifab_build_edge(chkdata_edge(n,i),mla%la(n),2,0,i) - call multifab_copy_c(chkdata_edge(n,i),1,umac(n,i),1,1) - call multifab_copy_c(chkdata_edge(n,i),2,umac_sum(n,i),1,1) + if (advection_type.eq.0) then + call multifab_build_edge(chkdata_edge(n,i),mla%la(n),(nspecies+1)*2+2,0,i) + call multifab_copy_c(chkdata_edge(n,i),1,umac(n,i),1,1) + call multifab_copy_c(chkdata_edge(n,i),2,umac_sum(n,i),1,1) + call multifab_copy_c(chkdata_edge(n,i),3,mass_fluxes(n,i),1,nspecies+1) + call multifab_copy_c(chkdata_edge(n,i),(nspecies+1)+3,mass_fluxes_sum(n,i),1,nspecies+1) + else + call multifab_build_edge(chkdata_edge(n,i),mla%la(n),2,0,i) + call multifab_copy_c(chkdata_edge(n,i),1,umac(n,i),1,1) + call multifab_copy_c(chkdata_edge(n,i),2,umac_sum(n,i),1,1) + endif endif end do end do + write(unit=check_index,fmt='(i8.8)') istep_to_write sd_name = trim(check_base_name) // check_index diff --git a/staggered_grid/src_lowMach/inhomogeneous_bc_val.f90 b/staggered_grid/src_lowMach/inhomogeneous_bc_val.f90 index ee64096a..8c4e6e19 100644 --- a/staggered_grid/src_lowMach/inhomogeneous_bc_val.f90 +++ b/staggered_grid/src_lowMach/inhomogeneous_bc_val.f90 @@ -8,13 +8,13 @@ module inhomogeneous_bc_val_module use probin_common_module, only: prob_lo, prob_hi, wallspeed_lo, wallspeed_hi, prob_type, & nspecies, algorithm_type, rho0, n_cells use probin_charged_module, only: Epot_wall, Epot_wall_bc_type, zero_charge_on_wall_type, bc_function_type, & - L_pos, L_trans, L_zero, induced_charge_eo, E_ext_value + L_pos, L_trans, L_zero, induced_charge_eo, E_ext_value, ac_iceo implicit none private - public :: scalar_bc, transport_bc, inhomogeneous_bc_val_2d, inhomogeneous_bc_val_3d + public :: scalar_bc, transport_bc, inhomogeneous_bc_val_2d, inhomogeneous_bc_val_3d, alternating_current_efield contains @@ -201,7 +201,11 @@ function inhomogeneous_bc_val_2d(comp,x,y,time_in) result(val) ! phi_periodic = x*E_0. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (induced_charge_eo) then - val = 1.d0*(x-(prob_lo(1)+prob_hi(1))/2.d0)*E_ext_value(1) + if (ac_iceo) then + val = 1.d0*(x-(prob_lo(1)+prob_hi(1))/2.d0)*alternating_current_efield(time) + else + val = 1.d0*(x-(prob_lo(1)+prob_hi(1))/2.d0)*E_ext_value(1) + endif else val = Epot_wall(1,2) end if @@ -218,41 +222,41 @@ function inhomogeneous_bc_val_2d(comp,x,y,time_in) result(val) else val = 0.d0 end if - else if (bc_function_type.eq.1) then ! Piecewise cubic potential--note we were lazy and assumed prob_lo(1) = 0.d0 - if (y .eq. prob_lo(2)) then - if (x .lt. L_pos) then ! we're in first positive region + else if (bc_function_type.eq.1) then ! Piecewise cubic potential--note we were lazy and assumed prob_lo(1) = 0.d0 + if (y .eq. prob_lo(2)) then ! lower boundary + if (x .lt. L_pos) then ! we're in first positive region val = Epot_wall(1,2) - else if (x .lt. (L_pos + L_trans)) then ! we're in first transition region + else if (x .lt. (L_pos + L_trans)) then ! we're in first transition region delta_x = prob_hi(1)/n_cells(1) s_j = x - L_pos val = 1.d0/delta_x*(0.5d0*Epot_wall(1,2)/L_trans**3*((s_j + delta_x/2.d0)**4 - (s_j - delta_x/2.d0)**4) - & Epot_wall(1,2)/L_trans**2*((s_j + delta_x/2.d0)**3 - (s_j - delta_x/2.d0)**3)) + Epot_wall(1,2) - else if (x .lt. (L_pos + L_trans + L_zero)) then ! we're in zero region + else if (x .lt. (L_pos + L_trans + L_zero)) then ! we're in zero region val = 0.d0 else if (x .lt. (L_pos + 2.d0*L_trans + L_zero)) then !we're in second transition region delta_x = prob_hi(1)/n_cells(1) s_j = x - (L_pos + L_trans + L_zero) val = 1.d0/delta_x*(-0.5d0*Epot_wall(1,2)/L_trans**3*((s_j + delta_x/2.d0)**4 - (s_j - delta_x/2.d0)**4) + & Epot_wall(1,2)/L_trans**2*((s_j + delta_x/2.d0)**3 - (s_j - delta_x/2.d0)**3)) - else ! we're in second positive region + else ! we're in second positive region val = Epot_wall(1,2) end if - else if (y .eq. prob_hi(2)) then - if (x .lt. L_pos) then ! we're in first positive region + else if (y .eq. prob_hi(2)) then ! upper boundary + if (x .lt. L_pos) then ! we're in first positive region val = Epot_wall(2,2) - else if (x .lt. (L_pos + L_trans)) then ! we're in first transition region + else if (x .lt. (L_pos + L_trans)) then ! we're in first transition region delta_x = prob_hi(1)/n_cells(1) s_j = x - L_pos val = 1.d0/delta_x*(0.5d0*Epot_wall(2,2)/L_trans**3*((s_j + delta_x/2.d0)**4 - (s_j - delta_x/2.d0)**4) - & Epot_wall(2,2)/L_trans**2*((s_j + delta_x/2.d0)**3 - (s_j - delta_x/2.d0)**3)) + Epot_wall(2,2) - else if (x .lt. (L_pos + L_trans + L_zero)) then ! we're in zero region + else if (x .lt. (L_pos + L_trans + L_zero)) then ! we're in zero region val = 0.d0 else if (x .lt. (L_pos + 2.d0*L_trans + L_zero)) then !we're in second transition region delta_x = prob_hi(1)/n_cells(1) s_j = x - (L_pos + L_trans + L_zero) val = 1.d0/delta_x*(-0.5d0*Epot_wall(2,2)/L_trans**3*((s_j + delta_x/2.d0)**4 - (s_j - delta_x/2.d0)**4) + & Epot_wall(2,2)/L_trans**2*((s_j + delta_x/2.d0)**3 - (s_j - delta_x/2.d0)**3)) - else ! we're in second positive region + else ! we're in second positive region val = Epot_wall(2,2) end if else @@ -271,6 +275,39 @@ function inhomogeneous_bc_val_2d(comp,x,y,time_in) result(val) end function inhomogeneous_bc_val_2d + ! This function allows us to specify an external electric field that oscillates in time. + ! The current field is a mollified square wave made by stitching together hyperbolic tangents. + ! + ! Currently the period and frequency of the wave is hard-coded in this function. + function alternating_current_efield(time) result(val) + + ! inputs/output + real(kind=dp_t), intent(in) :: time + real(kind=dp_t) :: val + + ! local variables + real(kind=dp_t) :: omega, T, shift ! omega controls how quickly the applied field goes from + ! positive to negative, and vice-versa. The larger the value, + ! the less quickly it transitions. + ! T is the period of oscillation. + integer :: branch + + ! initialize T and omega + T = 1.0d-4 ! these values were selected for a proof of concept AC-ICEO simulation. + omega = 0.05d0*T + + + branch = modulo(int(floor(time/T)), 2) !tells us if we use tanh or -1*tanh + shift = floor(time/T)*T + + if (branch.eq.0) then + val = -1.d0*E_ext_value(1)*tanh(((time-shift) - T/2.d0)/omega) + else + val = E_ext_value(1)*tanh(((time-shift) - T/2.d0)/omega) + endif + + end function alternating_current_efield + function inhomogeneous_bc_val_3d(comp,x,y,z,time_in) result(val) integer , intent(in ) :: comp diff --git a/staggered_grid/src_lowMach/init_lowmach.f90 b/staggered_grid/src_lowMach/init_lowmach.f90 index c3e7d8d6..a92636d5 100644 --- a/staggered_grid/src_lowMach/init_lowmach.f90 +++ b/staggered_grid/src_lowMach/init_lowmach.f90 @@ -208,7 +208,7 @@ subroutine init_rho_and_umac(mla,rho,umac,dx,time,the_bc_level) ! c1,c2,c3,u if (prob_type.eq.17) then if (nspecies.ne.3) then - call bl_error("Only nspecies = 3 is currently supported for prob_type 17") + call bl_error("Only nspecies = 3 is currently supported for prob_type=17") endif open(unit=34,file="c1_1d_vals.txt") open(unit=35,file="c2_1d_vals.txt") @@ -224,7 +224,12 @@ subroutine init_rho_and_umac(mla,rho,umac,dx,time,the_bc_level) close(35) close(36) close(37) - end if + else + c1_1d_arr=0.0d0 + c2_1d_arr=0.0d0 + c3_1d_arr=0.0d0 + u_1d_arr=0.0d0 + end if ! looping over boxes do n=1,nlevs @@ -647,7 +652,7 @@ subroutine init_rho_and_umac_2d(c,ng_c,rho,ng_r,u,v,ng_u,lo,hi,dx,time, & enddo enddo - case (15,16) + case (15) !========================================================= ! Discontinuous band in central 1/2 (type=15) or 1/3 (type=16) of domain diff --git a/staggered_grid/src_lowMach/main_driver.f90 b/staggered_grid/src_lowMach/main_driver.f90 index aa3fc660..bea9c39f 100644 --- a/staggered_grid/src_lowMach/main_driver.f90 +++ b/staggered_grid/src_lowMach/main_driver.f90 @@ -40,15 +40,17 @@ subroutine main_driver() use user_analysis use probin_common_module, only: prob_lo, prob_hi, n_cells, dim_in, hydro_grid_int, & max_grid_size, n_steps_save_stats, n_steps_skip, & + reset_averages, & plot_int, chk_int, seed, stats_int, bc_lo, bc_hi, restart, & probin_common_init, print_int, nspecies, dx_saved, & advection_type, fixed_dt, dt_saved, max_step, cfl, & algorithm_type, variance_coef_mom, initial_variance_mom, & variance_coef_mass, barodiffusion_type, use_bl_rng, & density_weights, rhobar, rho0, analyze_conserved, rho_eos_form, & - molmass, k_B, reset_tavg_vals, reset_tavg_step, stat_save_type, analyze_cuts + molmass, k_B, stat_save_type, analyze_cuts - use probin_multispecies_module, only: Dbar, start_time, probin_multispecies_init, c_init, T_init, c_bc + use probin_multispecies_module, only: Dbar, start_time, probin_multispecies_init, c_init, & + T_init, c_bc, additive_noise use probin_gmres_module, only: probin_gmres_init @@ -66,7 +68,7 @@ subroutine main_driver() ! quantities will be allocated with (nlevs,dm) components real(kind=dp_t), allocatable :: dx(:,:) real(kind=dp_t) :: dt,time,runtime1,runtime2,Dbar_max,dt_diffusive - integer :: n,nlevs,i,dm,istep,ng_s,init_step,n_Dbar,m + integer :: n,nlevs,i,dm,istep,ng_s,init_step,n_Dbar,m,comp type(box) :: bx type(ml_boxarray) :: mba type(ml_layout) :: mla @@ -106,6 +108,10 @@ subroutine main_driver() type(multifab), allocatable :: chem_rate(:) + ! these are only for the inertial algorithm + type(multifab), allocatable :: mass_fluxes(:,:) + type(multifab), allocatable :: charge_fluxes(:,:) + ! for writing time-averaged umac to plotfile type(multifab), allocatable :: umac_sum(:,:) type(multifab), allocatable :: umac_avg(:,:) @@ -118,6 +124,12 @@ subroutine main_driver() type(multifab), allocatable :: Epot_sum(:) type(multifab), allocatable :: Epot_avg(:) + ! for writing time-averaged charge fluxes to plotfile, only if inertial algorithm is used + type(multifab), allocatable :: mass_fluxes_sum(:,:) + type(multifab), allocatable :: mass_fluxes_avg(:,:) + type(multifab), allocatable :: charge_fluxes_sum(:,:) + type(multifab), allocatable :: charge_fluxes_avg(:,:) + ! for algorithm_type=6, return total mass fluxes (diff + stoch + adv) type(multifab), allocatable :: total_mass_flux(:,:) @@ -129,6 +141,7 @@ subroutine main_driver() ! misc real(kind=dp_t) :: max_charge, max_charge_abs, debye_len, sum_of_boundary_lens, delta_x + logical :: compute_and_save ! DONEV FIXME real(kind=dp_t) :: rho_temp, w_temp(1:5), w_mol(1:3) @@ -136,6 +149,8 @@ subroutine main_driver() m_Cl=5.887108600000000d-023, & m_H =1.673723600000000d-024, & m_OH=2.824068560000000d-023 + + integer :: num_avg_snapshots !============================================================== ! Initialization @@ -207,6 +222,12 @@ subroutine main_driver() allocate(Temp(nlevs)) allocate(diff_mass_fluxdiv(nlevs),stoch_mass_fluxdiv(nlevs)) allocate(stoch_mass_flux(nlevs,dm)) + allocate(mass_fluxes(nlevs,dm)) + allocate(mass_fluxes_sum(nlevs,dm)) + allocate(mass_fluxes_avg(nlevs,dm)) + allocate(charge_fluxes(nlevs,dm)) + allocate(charge_fluxes_sum(nlevs,dm)) + allocate(charge_fluxes_avg(nlevs,dm)) allocate(umac(nlevs,dm),mtemp(nlevs,dm),rhotot_fc(nlevs,dm),gradp_baro(nlevs,dm)) allocate(umac_sum(nlevs,dm)) allocate(umac_avg(nlevs,dm)) @@ -252,18 +273,88 @@ subroutine main_driver() init_step = restart + 1 - ! build the ml_layout ! read in time and dt from checkpoint - ! build and fill rho, rhotot, pi, umac, and umac_sum + ! build and fill rho, rhotot, pi, umac, umac_sum, + ! and mass_fluxes and mass_fluxes_sum if advection_type = 0 + ! and Epot, Epot_sum, and grad_Epot_old if use_charged_fluid is on call initialize_from_restart(mla,time,dt,rho_old,rho_sum,rhotot_old,pi,umac,umac_sum,pmask, & - Epot,Epot_sum,grad_Epot_old) + mass_fluxes,mass_fluxes_sum,Epot,Epot_sum,grad_Epot_old) ! enabled a fixed_dt that is different from the previous run if (fixed_dt .gt. 0.d0) then dt = fixed_dt end if + + ! upon restart, if n_steps_skip is greater than or equal to the restart number, zero the + ! sums used to compute time-averaged quantities + ! this may be redundent if n_steps_skip was not reached in the earlier simulation + ! (since the sum are already zero) but it allows us to reset the sum on restart if + ! we retroactively decide to increase n_steps_skip + if (n_steps_skip .ge. restart) then + ! reset rho_sum, Epot_sum, umac_sum, and mass_fluxes_sum + do n=1,nlevs + call setval(rho_sum(n),0.d0) + if (use_charged_fluid) then + call setval(Epot_sum(n),0.d0, all=.true.) + endif + do i=1,dm + call setval(umac_sum(n,i),0.d0) + end do + if (advection_type .eq. 0) then + do i=1,dm + call setval(mass_fluxes_sum(n,i),0.d0) + end do + end if + end do + end if + + ! Now we have mass_fluxes_sum (if advection_type = 0) and if use_charged_fluid, + ! then we'll build charge_fluxes and charge_fluxes_sum + ! and then convert mass_fluxes to charge_fluxes + ! + ! First copy mass fluxes into charge fluxes multifab. + ! Then convert the mass fluxes to charge fluxes. + if (use_charged_fluid .and. algorithm_type.eq.0) then + + ! build charge fluxes + do n=1,nlevs + do i=1,dm + call multifab_build_edge(charge_fluxes(n,i),mla%la(n),nspecies+1,0,i) + call multifab_build_edge(charge_fluxes_sum(n,i),mla%la(n),nspecies+1,0,i) + enddo + enddo + ! populate charge_fluxes w/ mass fluxes + do n=1,nlevs + do i=1,dm + call multifab_copy_c(charge_fluxes(n,i),1,mass_fluxes(n,i),1,nspecies,0) + call multifab_copy_c(charge_fluxes_sum(n,i),1,mass_fluxes_sum(n,i),1,nspecies,0) + enddo + enddo + + ! convert each mass flux to a charge flux + do n=1,nlevs + do i=1,dm + do comp=1,nspecies + call multifab_mult_mult_s_c(charge_fluxes(n,i),comp,charge_per_mass(comp),1,0) + call multifab_mult_mult_s_c(charge_fluxes_sum(n,i),comp,charge_per_mass(comp),1,0) + enddo + enddo + enddo + + ! now add up each individual charge flux for the total charge flux + do n=1,nlevs + do i=1,dm + call multifab_copy_c(charge_fluxes(n,i),nspecies+1,charge_fluxes(n,i),1,1,0) + call multifab_copy_c(charge_fluxes_sum(n,i),nspecies+1,charge_fluxes_sum(n,i),1,1,0) + do comp=2,nspecies + call multifab_plus_plus_c(charge_fluxes(n,i),nspecies+1,charge_fluxes(n,i),comp,1,0) + call multifab_plus_plus_c(charge_fluxes_sum(n,i),nspecies+1,charge_fluxes_sum(n,i),comp,1,0) + enddo + enddo + enddo + endif else init_step = 1 @@ -319,7 +410,6 @@ subroutine main_driver() call setval(umac_sum(n,i),0.d0) end do end do - end if ! data structures to help with reservoirs @@ -417,28 +507,28 @@ subroutine main_driver() dx_in=dx(n,:)) end do -if(.false.) then - ! DONEV FIXME temporary debugging: - w_mol=(/0.0014d0,0.0053d0,0.0037d0/) ! Random values - if(nspecies==5) then - ! 1=HCl, 2=NaOH, 3=NaCl, 4=H2O - ! 1=Na+, 2=Cl-, 3=H+, 4=OH-, 5=H2O - ! w_Na = m_Na / (m_Na+m_OH) * w_NaOH + m_Na / (m_Na+m_Cl) * w_NaCl - w_temp(1) = m_Na / (m_Na+m_OH) * w_mol(2) + m_Na / (m_Na+m_Cl) * w_mol(3) - ! w_Cl = m_Cl / (m_H+m_Cl) * w_HCl + m_Cl / (m_Na+m_Cl) * w_NaCl - w_temp(2) = m_Cl / (m_H+m_Cl) * w_mol(1) + m_Cl / (m_Na+m_Cl) * w_mol(3) - ! w_H = m_H / (m_H+m_Cl) * w_HCl - w_temp(3) = m_H / (m_H+m_Cl) * w_mol(1) - ! w_OH = m_OH / (m_Na+m_OH) * w_NaOH - w_temp(4) = m_OH / (m_Na+m_OH) * w_mol(2) - else - w_temp(1:3)=w_mol - end if - w_temp(nspecies) = 1.0d0 - sum(w_temp(1:nspecies)) - call compute_rho_eos(w_temp(1:nspecies)*rho0, rho_temp) - write(*,*) " w_H2O=", w_temp(nspecies), " rho_eos=", rho_temp, " w=", w_temp(1:nspecies-1) - stop -end if + if(.false.) then + ! DONEV FIXME temporary debugging: + w_mol=(/0.0014d0,0.0053d0,0.0037d0/) ! Random values + if(nspecies==5) then + ! 1=HCl, 2=NaOH, 3=NaCl, 4=H2O + ! 1=Na+, 2=Cl-, 3=H+, 4=OH-, 5=H2O + ! w_Na = m_Na / (m_Na+m_OH) * w_NaOH + m_Na / (m_Na+m_Cl) * w_NaCl + w_temp(1) = m_Na / (m_Na+m_OH) * w_mol(2) + m_Na / (m_Na+m_Cl) * w_mol(3) + ! w_Cl = m_Cl / (m_H+m_Cl) * w_HCl + m_Cl / (m_Na+m_Cl) * w_NaCl + w_temp(2) = m_Cl / (m_H+m_Cl) * w_mol(1) + m_Cl / (m_Na+m_Cl) * w_mol(3) + ! w_H = m_H / (m_H+m_Cl) * w_HCl + w_temp(3) = m_H / (m_H+m_Cl) * w_mol(1) + ! w_OH = m_OH / (m_Na+m_OH) * w_NaOH + w_temp(4) = m_OH / (m_Na+m_OH) * w_mol(2) + else + w_temp(1:3)=w_mol + end if + w_temp(nspecies) = 1.0d0 - sum(w_temp(1:nspecies)) + call compute_rho_eos(w_temp(1:nspecies)*rho0, rho_temp) + write(*,*) " w_H2O=", w_temp(nspecies), " rho_eos=", rho_temp, " w=", w_temp(1:nspecies-1) + stop + end if !======================================================= ! Build multifabs for all the variables @@ -494,6 +584,43 @@ subroutine main_driver() end do end if + ! Only build the mass flux to be kept track of below if advection type is centered explicit + ! + ! mass_fluxes contains nspecies+1 components--one for each species plus a component for + ! all the species summed up. + ! + ! If use_charged_fluid is on, then we'll keep track of charge_fluxes as well, + ! where charge_flux = mass_flux * (charge/mass) + if (advection_type.eq.0) then + do n=1,nlevs + do i=1,dm + + ! always build the averages + call multifab_build_edge(mass_fluxes_avg(n,i),mla%la(n),nspecies+1,0,i) + call setval(mass_fluxes_avg(n,i),0.d0) + if (use_charged_fluid) then + call multifab_build_edge(charge_fluxes_avg(n,i),mla%la(n),nspecies+1,0,i) + call setval(charge_fluxes_avg(n,i),0.d0) + endif + + ! only build the mass/charge fluxes and their sums if a checkpoint file wasn't provided + if (restart .lt. 0) then + call multifab_build_edge(mass_fluxes(n,i),mla%la(n),nspecies+1,1,i) + call multifab_build_edge(mass_fluxes_sum(n,i),mla%la(n),nspecies+1,0,i) + call setval(mass_fluxes(n,i),0.d0) + call setval(mass_fluxes_sum(n,i),0.d0) + if (use_charged_fluid) then + call multifab_build_edge(charge_fluxes(n,i),mla%la(n),nspecies+1,1,i) + call multifab_build_edge(charge_fluxes_sum(n,i),mla%la(n),nspecies+1,0,i) + call setval(charge_fluxes(n,i),0.d0) + call setval(charge_fluxes_sum(n,i),0.d0) + end if + endif + + end do + end do + end if + if (use_charged_fluid) then do n=1,nlevs call multifab_build(charge_old(n) ,mla%la(n),1,1) @@ -553,6 +680,20 @@ subroutine main_driver() ! Initialize values !===================================================================== + ! compute sqrtLonsgaer for additive noise + if (additive_noise) then + ! use rho_new as a temporary holding spot for the 'equilibrium' condition + call init_additive_sqrtLonsager(mla) + ! initialize rho and umac in valid region only + call init_rho_and_umac(mla,rho_new,umac,dx,time,the_bc_tower%bc_tower_array) + ! compute rhotot from rho in VALID REGION + call compute_rhotot(mla,rho_new,rhotot_new) + ! fill rho and rhotot ghost cells + call fill_rho_rhotot_ghost(mla,rho_new,rhotot_new,dx,the_bc_tower) + ! computer noise multiplier + call compute_sqrtLonsager_fc(mla,rho_new,rhotot_new,additive_sqrtLonsager_fc,dx) + end if + if (use_charged_fluid) then ! set these to zero @@ -576,6 +717,7 @@ subroutine main_driver() end do end do + if(electroneutral) then call dot_with_z(mla,rho_old,charge_old,abs_z=.true.) max_charge_abs = multifab_norm_inf(charge_old(1)) ! This will be saved for reuse @@ -587,9 +729,13 @@ subroutine main_driver() call dot_with_z(mla,rho_old,charge_old) ! multiply by total volume (all 3 dimensions, even for 2D problems) - if (use_charged_fluid) then ! NOTE: we are using rho = 1 here, so the below is a close approximation to debye length - debye_len =sqrt(dielectric_const*k_B*T_init(1)/ & - (rho0*sum(c_init(1,1:nspecies)*molmass(1:nspecies)*charge_per_mass(1:nspecies)**2))) + if (use_charged_fluid.and.print_debye_len) then ! NOTE: we are using rho = 1 here, so the below is a close approximation to debye len + debye_len = 0.d0 + do comp=1,nspecies + debye_len = debye_len + c_init(1,comp)*molmass(comp)*charge_per_mass(comp)*charge_per_mass(comp) + enddo + debye_len = sqrt(dielectric_const*k_B*T_init(1)/debye_len) + if (parallel_IOprocessor()) then print*, 'Debye length $\lambda_D$ is approx: ', debye_len endif @@ -749,13 +895,26 @@ subroutine main_driver() end if ! We will also pass temperature - call initialize_hydro_grid(mla,rho_old,dt,dx,namelist_file=un, & - nspecies_in=nspecies, & - nscal_in=nscal, & - exclude_last_species_in=.false., & - analyze_velocity=.true., & - analyze_density=.true., & - analyze_temperature=.true.) + if (use_charged_fluid) then + + call initialize_hydro_grid(mla,rho_old,dt,dx,namelist_file=un, & + nspecies_in=nspecies, & + nscal_in=nscal, & + exclude_last_species_in=.false., & + analyze_velocity=.true., & + analyze_density=.true., & + analyze_temperature=.true., & + analyze_charge_fluxes=.true.) + else + call initialize_hydro_grid(mla,rho_old,dt,dx,namelist_file=un, & + nspecies_in=nspecies, & + nscal_in=nscal, & + exclude_last_species_in=.false., & + analyze_velocity=.true., & + analyze_density=.true., & + analyze_temperature=.true., & + analyze_charge_fluxes=.false.) + endif close(unit=un) @@ -830,19 +989,27 @@ subroutine main_driver() ! Hydrogrid analysis and output for initial data !===================================================================== - do n=1,nlevs - do i=1,dm - call multifab_copy_c(umac_avg(n,i),1,umac(n,i),1,1,0) + if (n_steps_skip .eq. 0) then + do n=1,nlevs + do i=1,dm + call multifab_copy_c(umac_avg(n,i),1,umac(n,i),1,1,0) + end do end do - end do + else + do n=1,nlevs + do i=1,dm + call multifab_setval(umac_avg(n,i),0.d0,all=.true.) + end do + end do + end if ! write initial plotfile if (plot_int .gt. 0) then if (parallel_IOProcessor()) then write(*,*) 'writing initial plotfile 0' end if - call write_plotfile(mla,rho_old,rho_avg,rhotot_old,Temp,umac,umac_avg,pi,Epot, & - Epot_avg,grad_Epot_old,gradPhiApprox,0,dx,time) + call write_plotfile(mla,rho_old,rho_avg,rhotot_old,Temp,umac,umac_avg,pi,mass_fluxes,mass_fluxes_avg,Epot, & + Epot_avg,grad_Epot_old,gradPhiApprox,charge_fluxes,charge_fluxes_avg,0,dx,time) end if ! write initial checkpoint @@ -850,13 +1017,14 @@ subroutine main_driver() if (parallel_IOProcessor()) then write(*,*) 'writing initial checkpoint 0' end if - call checkpoint_write(mla,rho_old,rho_sum,rhotot_old,pi,umac,umac_sum,Epot,Epot_sum,grad_Epot_old,time,dt,0) + call checkpoint_write(mla,rho_old,rho_sum,rhotot_old,pi,umac,umac_sum,mass_fluxes,mass_fluxes_sum, & + Epot,Epot_sum,grad_Epot_old,time,dt,0) end if if (stats_int .gt. 0) then ! write initial vertical and horizontal averages (hstat and vstat files) if (use_charged_fluid) then - call print_stats(mla,dx,0,time,umac=umac,rho=rho_old,temperature=Temp,scalars=Epot) + call print_stats(mla,dx,0,time,umac=umac,rho=rho_old,temperature=Temp,scalars=Epot,charge_fluxes=charge_fluxes) else call print_stats(mla,dx,0,time,umac=umac,rho=rho_old,temperature=Temp) end if @@ -867,7 +1035,12 @@ subroutine main_driver() ! Add this snapshot to the average in HydroGrid if (hydro_grid_int > 0) then - call analyze_hydro_grid(mla,dt,dx,istep,umac=umac,rho=rho_old,temperature=Temp) + if (use_charged_fluid) then + call analyze_hydro_grid(mla,dt,dx,0,umac=umac,rho=rho_old,temperature=Temp, & + scalars=Epot) + else + call analyze_hydro_grid(mla,dt,dx,0,umac=umac,rho=rho_old,temperature=Temp) + end if end if if (hydro_grid_int > 0 .and. n_steps_save_stats > 0) then @@ -887,7 +1060,7 @@ subroutine main_driver() ! Begin time stepping loop !======================================================= - do istep=init_step,max_step + TimeLoop: do istep=init_step,max_step if (fixed_dt .le. 0.d0) then call estdt(mla,umac,dx,dt) @@ -898,8 +1071,10 @@ subroutine main_driver() end if if (parallel_IOProcessor()) then - if ( (print_int .gt. 0 .and. mod(istep,print_int) .eq. 0) ) & - print*,"Begin Advance; istep =",istep,"dt =",dt,"time =",time + if ( print_int .gt. 0 ) then + if (mod(istep,print_int) .eq. 0) & + print*,"Begin Advance; istep =",istep,"dt =",dt,"time =",time + end if end if runtime1 = parallel_wtime() @@ -913,14 +1088,18 @@ subroutine main_driver() ! diff/stoch_mass_fluxdiv could be built locally within the overdamped ! routine, but since we have them around anyway for inertial we pass them in if (algorithm_type .eq. 0) then - ! algorithm_type=0: inertial + ! algorithm_type=0: inertial !SC + ! in this case, we record the species' mass fluxes. After this function call, + ! mass_fluxes is actually each individual species call advance_timestep_inertial(mla,umac,rho_old,rho_new,rhotot_old,rhotot_new, & gradp_baro,pi,eta,eta_ed,kappa,Temp,Temp_ed, & diff_mass_fluxdiv, stoch_mass_fluxdiv,stoch_mass_flux, & + mass_fluxes, & dx,dt,time,the_bc_tower,istep, & grad_Epot_old,grad_Epot_new, & charge_old,charge_new,Epot, & permittivity) + else if (algorithm_type .eq. 2) then ! algorithm_type=2: overdamped with 2 RNG call advance_timestep_overdamped(mla,umac,rho_old,rho_new,rhotot_old,rhotot_new, & @@ -991,88 +1170,143 @@ subroutine main_driver() print*,'' end if end if - end if - - ! for writing time-averaged umac, rho, and Epot to plotfile - if (istep .gt. n_steps_skip) then - ! Note: reset time avg quantities is only possible if reset_tavg_step >= n_steps_skip - ! Also note: if reset is turned ON, and istep is between n_steps_skip and reset_tavg_step, - ! then the code below will not track the average (since it will get reset anyways) - ! - ! TL;DR: don't turn reset_tavg_vals ON unless you want to use it. - if (reset_tavg_vals .and. (reset_tavg_step.ge.n_steps_skip)) then - if (istep.eq.reset_tavg_step) then - do n=1,nlevs - ! reset rho_sum, epot_sum, and umac_sum to be 0 - call setval(rho_sum(n),0.d0) - call setval(Epot_sum(n),0.d0, all=.true.) - do i=1,dm - call setval(umac_sum(n,i),0.d0) - end do - end do - else if (istep .gt. reset_tavg_step) then - ! do the normal averaging, starting from the reset step - do n=1,nlevs - ! first do rho - call multifab_plus_plus_c(rho_sum(n),1,rho_new(n),1,nspecies,0) - call multifab_copy_c(rho_avg(n),1,rho_sum(n),1,nspecies,0) - call multifab_mult_mult_s_c(rho_avg(n),1,(1.d0/(istep-reset_tavg_step)),nspecies,0) - - ! next do Epot - if (use_charged_fluid) then - call multifab_plus_plus_c(Epot_sum(n),1,Epot(n),1,1,0) - call multifab_copy_c(Epot_avg(n),1,Epot_sum(n),1,1,0) - call multifab_mult_mult_s_c(Epot_avg(n),1,(1.d0/(istep-reset_tavg_step)),1,0) - end if - - ! lastly do umac - do i=1,dm - call multifab_plus_plus_c(umac_sum(n,i),1,umac(n,i),1,1,0) - call multifab_copy_c(umac_avg(n,i),1,umac_sum(n,i),1,1,0) - call multifab_mult_mult_s_c(umac_avg(n,i),1,(1.d0/(istep-reset_tavg_step)),1,0) - end do - end do - end if - else - do n=1,nlevs - ! first do rho - call multifab_plus_plus_c(rho_sum(n),1,rho_new(n),1,nspecies,0) - call multifab_copy_c(rho_avg(n),1,rho_sum(n),1,nspecies,0) - call multifab_mult_mult_s_c(rho_avg(n),1,(1.d0/(istep-n_steps_skip)),nspecies,0) + end if + + ! If use charged fluid and inertial algorithm, first copy mass fluxes into charge fluxes multifab. + ! Then convert the mass fluxes to charge fluxes. + if (use_charged_fluid.and.(algorithm_type.eq.0)) then + ! populate charge_fluxes w/ mass fluxes + do n=1,nlevs + do i=1,dm + call multifab_copy_c(charge_fluxes(n,i),1,mass_fluxes(n,i),1,nspecies,0) + enddo + enddo + + ! convert each mass flux to a charge flux + do n=1,nlevs + do i=1,dm + do comp=1,nspecies + call multifab_mult_mult_s_c(charge_fluxes(n,i),comp,charge_per_mass(comp),1,0) + enddo + enddo + enddo + + ! now add up each individual charge flux for the total charge flux + do n=1,nlevs + do i=1,dm + call multifab_copy_c(charge_fluxes(n,i),nspecies+1,charge_fluxes(n,i),1,1,0) + do comp=2,nspecies + call multifab_plus_plus_c(charge_fluxes(n,i),nspecies+1,charge_fluxes(n,i),comp,1,0) + enddo + enddo + enddo + endif - ! next do Epot + ! replace .true. with new reset_averages flag + if (reset_averages) then + + ! reset averages every n_steps_skip+1 + if(n_steps_skip>0) then + if (mod(istep,n_steps_skip) .eq. 1) then + do n=1,nlevs + call setval(rho_sum(n),0.d0) if (use_charged_fluid) then - call multifab_plus_plus_c(Epot_sum(n),1,Epot(n),1,1,0) - call multifab_copy_c(Epot_avg(n),1,Epot_sum(n),1,1,0) - call multifab_mult_mult_s_c(Epot_avg(n),1,(1.d0/(istep-n_steps_skip)),1,0) + call setval(Epot_sum(n),0.d0) end if - - ! lastly do umac do i=1,dm - call multifab_plus_plus_c(umac_sum(n,i),1,umac(n,i),1,1,0) - call multifab_copy_c(umac_avg(n,i),1,umac_sum(n,i),1,1,0) - call multifab_mult_mult_s_c(umac_avg(n,i),1,(1.d0/(istep-n_steps_skip)),1,0) + call setval(umac_sum(n,i),0.d0) + if (advection_type.eq.0) then + call setval(mass_fluxes_sum(n,i),0.d0) + if (use_charged_fluid) then + call setval(charge_fluxes_sum(n,i),0.d0) + end if + end if end do end do end if + end if + + if(n_steps_skip>0) then + if (mod(istep,n_steps_skip) .eq. 0) then + num_avg_snapshots = n_steps_skip + else + num_avg_snapshots = mod(istep,n_steps_skip) + end if + end if + + else + + num_avg_snapshots = istep-n_steps_skip + end if + ! for writing time-averaged umac, rho, Epot, mass fluxes, and charge fluxes to plotfile + if (istep .gt. n_steps_skip .or. reset_averages) then + + do n=1,nlevs + ! first do rho + call multifab_plus_plus_c(rho_sum(n),1,rho_new(n),1,nspecies,0) + call multifab_copy_c(rho_avg(n),1,rho_sum(n),1,nspecies,0) + call multifab_mult_mult_s_c(rho_avg(n),1,(1.d0/num_avg_snapshots),nspecies,0) + + ! next do Epot + if (use_charged_fluid) then + call multifab_plus_plus_c(Epot_sum(n),1,Epot(n),1,1,0) + call multifab_copy_c(Epot_avg(n),1,Epot_sum(n),1,1,0) + call multifab_mult_mult_s_c(Epot_avg(n),1,(1.d0/num_avg_snapshots),1,0) + end if + + ! do umac + do i=1,dm + call multifab_plus_plus_c(umac_sum(n,i),1,umac(n,i),1,1,0) + call multifab_copy_c(umac_avg(n,i),1,umac_sum(n,i),1,1,0) + call multifab_mult_mult_s_c(umac_avg(n,i),1,(1.d0/num_avg_snapshots),1,0) + end do + + ! do mass fluxes + if (advection_type.eq.0) then + do i=1,dm + call multifab_plus_plus_c(mass_fluxes_sum(n,i),1,mass_fluxes(n,i) ,1,nspecies+1,0) + call multifab_copy_c( mass_fluxes_avg(n,i),1,mass_fluxes_sum(n,i),1,nspecies+1,0) + call multifab_mult_mult_s_c(mass_fluxes_avg(n,i),1,(1.d0/num_avg_snapshots),nspecies+1,0) + end do + + ! do charge fluxes + if (use_charged_fluid) then + do i=1,dm + call multifab_plus_plus_c(charge_fluxes_sum(n,i),1,charge_fluxes(n,i) ,1,nspecies+1,0) + call multifab_copy_c( charge_fluxes_avg(n,i),1,charge_fluxes_sum(n,i),1,nspecies+1,0) + call multifab_mult_mult_s_c(charge_fluxes_avg(n,i),1,(1.d0/num_avg_snapshots),nspecies+1,0) + enddo + endif + endif + end do + + end if + + time = time + dt - if (parallel_IOProcessor()) then - if ( (print_int .gt. 0 .and. mod(istep,print_int) .eq. 0) ) & - print*,"End Advance; istep =",istep,"DT =",dt,"TIME =",time + if (parallel_IOProcessor() .and. (print_int .gt. 0) ) then + if ( mod(istep,print_int) .eq. 0) & + print*,"End Advance; istep =",istep,"DT =",dt,"TIME =",time end if runtime2 = parallel_wtime() - runtime1 call parallel_reduce(runtime1, runtime2, MPI_MAX, proc=parallel_IOProcessorNode()) - if (parallel_IOProcessor()) then - if ( (print_int .gt. 0 .and. mod(istep,print_int) .eq. 0) ) & + if (parallel_IOProcessor() .and. (print_int .gt. 0) ) then + if ( mod(istep,print_int) .eq. 0 ) & print*,'Time to advance timestep: ',runtime1,' seconds' end if - if ( (print_int .gt. 0 .and. mod(istep,print_int) .eq. 0) .or. & - (istep .eq. max_step) ) then + ! Compute conserved quantities periodically: + compute_and_save = .false. + if (istep .eq. max_step) compute_and_save = .true. + if ( (print_int .gt. 0) ) then + if ( (mod(istep,print_int) .eq. 0) ) compute_and_save = .true. + end if + if(compute_and_save) then + if (parallel_IOProcessor()) write(*,*) "After time step ", istep, " t=", time call sum_mass(rho_new, istep) ! print out the total mass to check conservation ! compute rhotot on faces @@ -1100,64 +1334,84 @@ subroutine main_driver() endif end if end if - - end if ! write plotfile at specific intervals - if (plot_int.gt.0 .and. ( (mod(istep,plot_int).eq.0) .or. (istep.eq.max_step)) ) then + compute_and_save = .false. + if (istep .eq. max_step) compute_and_save = .true. + if ( (plot_int .gt. 0) ) then + if ( (mod(istep,plot_int) .eq. 0) ) compute_and_save = .true. + end if + if(compute_and_save) then + if (parallel_IOProcessor()) then write(*,*) 'writing plotfiles after timestep =', istep end if - call write_plotfile(mla,rho_new,rho_avg,rhotot_new,Temp,umac,umac_avg,pi,Epot, & - Epot_avg,grad_Epot_new,gradPhiApprox,istep,dx,time) + call write_plotfile(mla,rho_new,rho_avg,rhotot_new,Temp,umac,umac_avg,pi,mass_fluxes,mass_fluxes_avg,Epot, & + Epot_avg,grad_Epot_new,gradPhiApprox,charge_fluxes,charge_fluxes_avg,istep,dx,time) end if ! write checkpoint at specific intervals - if ((chk_int.gt.0 .and. mod(istep,chk_int).eq.0)) then + if ( chk_int.gt.0 ) then + if ( mod(istep,chk_int).eq.0 ) then if (parallel_IOProcessor()) then write(*,*) 'writing checkpoint after timestep =', istep end if - call checkpoint_write(mla,rho_new,rho_sum,rhotot_new,pi,umac,umac_sum,Epot,Epot_sum,grad_Epot_new,time,dt,istep) + call checkpoint_write(mla,rho_new,rho_sum,rhotot_new,pi,umac,umac_sum,mass_fluxes,mass_fluxes_sum, & + Epot,Epot_sum,grad_Epot_new,time,dt,istep) + end if end if ! print out projection (average) and variance - if ( (stats_int > 0) .and. & - (mod(istep,stats_int) .eq. 0) ) then + if ( stats_int > 0 ) then + if ( mod(istep,stats_int) .eq. 0 ) then if (stat_save_type.eq.0) then ! compute vertical/horizontal averages of instantaneous fields if (use_charged_fluid) then - call print_stats(mla,dx,istep,time,umac=umac,rho=rho_new,temperature=Temp,scalars=Epot) + call print_stats(mla,dx,istep,time,umac=umac,rho=rho_new,temperature=Temp, & + scalars=Epot,charge_fluxes=charge_fluxes) else call print_stats(mla,dx,istep,time,umac=umac,rho=rho_new,temperature=Temp) endif else ! compute vertical/horizontal averages of time averaged fields - if (istep.lt.n_steps_skip) then - call bl_error("If stat_save_type is \ne 0, you cannot call print_stats until n_steps_skip has passed.") - end if - if (use_charged_fluid) then - call print_stats(mla,dx,istep,time,umac=umac_avg,rho=rho_avg,temperature=Temp,scalars=Epot_avg) - else - call print_stats(mla,dx,istep,time,umac=umac_avg,rho=rho_avg,temperature=Temp) - endif + if (istep.le.n_steps_skip) then + if (parallel_IOProcessor()) then + write(*,*) "WARNING: stat_save_type != 0 and we have not surpassed n_steps_skip yet" + write(*,*) "Skipping the call to print_stats (to write hstat/vstat files)" + end if + else + if (use_charged_fluid) then + call print_stats(mla,dx,istep,time,umac=umac_avg,rho=rho_avg,temperature=Temp, & + scalars=Epot_avg,charge_fluxes=charge_fluxes_avg) + else + call print_stats(mla,dx,istep,time,umac=umac_avg,rho=rho_avg,temperature=Temp) + endif + end if end if end if + end if if (istep .ge. n_steps_skip) then ! Add this snapshot to the average in HydroGrid - if ( (hydro_grid_int > 0) .and. & - ( mod(istep,hydro_grid_int) .eq. 0 ) ) then - call analyze_hydro_grid(mla,dt,dx,istep,umac=umac,rho=rho_new,temperature=Temp) + if ( hydro_grid_int > 0 ) then + if ( mod(istep,hydro_grid_int) .eq. 0 ) then + if (use_charged_fluid) then + call analyze_hydro_grid(mla,dt,dx,istep,umac=umac,rho=rho_new,temperature=Temp, & + scalars=Epot) + else + call analyze_hydro_grid(mla,dt,dx,istep,umac=umac,rho=rho_new,temperature=Temp) + end if + end if end if - if ( (hydro_grid_int > 0) .and. & - (n_steps_save_stats > 0) .and. & - ( mod(istep,n_steps_save_stats) .eq. 0 ) ) then + if ( (hydro_grid_int > 0) .and. (n_steps_save_stats > 0) ) then + if ( mod(istep,n_steps_save_stats) .eq. 0 ) then call save_hydro_grid(id=istep/n_steps_save_stats, step=istep) end if + end if end if @@ -1177,7 +1431,7 @@ subroutine main_driver() end do end if - end do + end do TimeLoop !======================================================= ! Destroy multifabs and layouts @@ -1227,6 +1481,21 @@ subroutine main_driver() end do end if + if (advection_type.eq.0) then !SC + do n=1,nlevs + do i=1,dm + call multifab_destroy(mass_fluxes(n,i)) + call multifab_destroy(mass_fluxes_avg(n,i)) + call multifab_destroy(mass_fluxes_sum(n,i)) + if (use_charged_fluid) then + call multifab_destroy(charge_fluxes(n,i)) + call multifab_destroy(charge_fluxes_sum(n,i)) + call multifab_destroy(charge_fluxes_avg(n,i)) + endif + end do + end do + endif + if (use_charged_fluid) then do n=1,nlevs call multifab_destroy(charge_old(n)) @@ -1270,6 +1539,10 @@ subroutine main_driver() end do end do + if (additive_noise) then + call destroy_additive_sqrtLonsager(mla) + end if + deallocate(lo,hi,pmask) deallocate(rho_old,rhotot_old,pi) deallocate(rho_sum, rho_avg) @@ -1277,6 +1550,8 @@ subroutine main_driver() deallocate(Temp) deallocate(diff_mass_fluxdiv,stoch_mass_fluxdiv) deallocate(stoch_mass_flux) + deallocate(mass_fluxes) !SC + deallocate(charge_fluxes) !SC deallocate(umac,mtemp,rhotot_fc,gradp_baro) deallocate(eta,kappa) deallocate(eta_ed) diff --git a/staggered_grid/src_lowMach/restart.f90 b/staggered_grid/src_lowMach/restart.f90 index 7f5014be..51c241e7 100644 --- a/staggered_grid/src_lowMach/restart.f90 +++ b/staggered_grid/src_lowMach/restart.f90 @@ -22,7 +22,8 @@ module restart_module contains - subroutine initialize_from_restart(mla,time,dt,rho,rho_sum,rhotot,pi,umac,umac_sum,pmask,Epot,Epot_sum,grad_Epot) + subroutine initialize_from_restart(mla,time,dt,rho,rho_sum,rhotot,pi,umac,umac_sum,pmask,mass_fluxes,mass_fluxes_sum, & + Epot,Epot_sum,grad_Epot) type(ml_layout),intent(out) :: mla real(dp_t) , intent( out) :: time,dt @@ -34,6 +35,8 @@ subroutine initialize_from_restart(mla,time,dt,rho,rho_sum,rhotot,pi,umac,umac_s type(multifab), intent(inout) :: Epot_sum(:) type(multifab), intent(inout) :: umac(:,:) type(multifab), intent(inout) :: umac_sum(:,:) + type(multifab), intent(inout) :: mass_fluxes(:,:) + type(multifab), intent(inout) :: mass_fluxes_sum(:,:) type(multifab), intent(inout) :: grad_Epot(:,:) logical , intent(in ) :: pmask(:) @@ -85,6 +88,17 @@ subroutine initialize_from_restart(mla,time,dt,rho,rho_sum,rhotot,pi,umac,umac_s call setval(umac(n,i),0.d0,all=.true.) call setval(umac_sum(n,i),0.d0,all=.true.) + if (advection_type.eq.0) then + call multifab_build_edge(mass_fluxes(n,i), mla%la(n), nspecies+1, 1, i) + call multifab_build_edge(mass_fluxes_sum(n,i), mla%la(n), nspecies+1, 0, i) !SC: Q: what should num_ghost_cell be here?? + + ! with mixed boundary conditions some of the corner umac ghost cells that + ! never affect the solution aren't filled, causing segfaults on some compilers + ! this prevents segfaults + call setval(mass_fluxes(n,i),0.d0,all=.true.) + call setval(mass_fluxes_sum(n,i),0.d0,all=.true.) + endif + if (use_charged_fluid) then call multifab_build_edge(grad_Epot(n,i), mla%la(n), 1, 1, i) call setval(grad_Epot(n,i),0.d0,all=.true.) ! for ghost cells as well @@ -105,23 +119,73 @@ subroutine initialize_from_restart(mla,time,dt,rho,rho_sum,rhotot,pi,umac,umac_s endif end do + ! edge data + ! + ! if use_charged_fluid is on, we build grad_Epot + ! if advection_type = 0, we build mass fluxes (both instantaneous and sum) do n=1,nlevs - call multifab_copy_c(umac(n,1),1,chkdata_edgex(n),1,1) - call multifab_copy_c(umac(n,2),1,chkdata_edgey(n),1,1) - call multifab_copy_c(umac_sum(n,1),1,chkdata_edgex(n),2,1) - call multifab_copy_c(umac_sum(n,2),1,chkdata_edgey(n),2,1) - if (use_charged_fluid) then - call multifab_copy_c(grad_Epot(n,1),1,chkdata_edgex(n),3,1) - call multifab_copy_c(grad_Epot(n,2),1,chkdata_edgey(n),3,1) - endif - if (dm .eq. 3) then - call multifab_copy_c(umac(n,3),1,chkdata_edgez(n),1,1) - call multifab_copy_c(umac_sum(n,3),1,chkdata_edgez(n),2,1) - if (use_charged_fluid) then - call multifab_copy_c(grad_Epot(n,3),1,chkdata_edgez(n),3,1) + if (use_charged_fluid) then + if (advection_type.eq.0) then + call multifab_copy_c(umac(n,1), 1,chkdata_edgex(n),1, 1) + call multifab_copy_c(umac(n,2), 1,chkdata_edgey(n),1, 1) + call multifab_copy_c(umac_sum(n,1), 1,chkdata_edgex(n),2, 1) + call multifab_copy_c(umac_sum(n,2), 1,chkdata_edgey(n),2, 1) + + call multifab_copy_c(mass_fluxes(n,1), 1,chkdata_edgex(n),3, nspecies+1) + call multifab_copy_c(mass_fluxes(n,2), 1,chkdata_edgey(n),3, nspecies+1) + call multifab_copy_c(mass_fluxes_sum(n,1),1,chkdata_edgex(n),(nspecies+1)+3, nspecies+1) + call multifab_copy_c(mass_fluxes_sum(n,2),1,chkdata_edgey(n),(nspecies+1)+3, nspecies+1) + call multifab_copy_c(grad_Epot(n,1), 1,chkdata_edgex(n),(nspecies+1)*2+3,1) + call multifab_copy_c(grad_Epot(n,2), 1,chkdata_edgey(n),(nspecies+1)*2+3,1) + if (dm .eq. 3) then + call multifab_copy_c(umac(n,3), 1,chkdata_edgez(n),1, 1) + call multifab_copy_c(umac_sum(n,3), 1,chkdata_edgez(n),2, 1) + call multifab_copy_c(mass_fluxes(n,3), 1,chkdata_edgez(n),3, nspecies+1) + call multifab_copy_c(mass_fluxes_sum(n,3),1,chkdata_edgez(n),(nspecies+1)+3, nspecies+1) + call multifab_copy_c(grad_Epot(n,3), 1,chkdata_edgez(n),(nspecies+1)*2+3,1) + end if + else ! advection_type \ne 0 + call multifab_copy_c(umac(n,1), 1,chkdata_edgex(n),1, 1) + call multifab_copy_c(umac(n,2), 1,chkdata_edgey(n),1, 1) + call multifab_copy_c(umac_sum(n,1), 1,chkdata_edgex(n),2, 1) + call multifab_copy_c(umac_sum(n,2), 1,chkdata_edgey(n),2, 1) + call multifab_copy_c(grad_Epot(n,1), 1,chkdata_edgex(n),3, 1) + call multifab_copy_c(grad_Epot(n,2), 1,chkdata_edgey(n),3, 1) + if (dm .eq. 3) then + call multifab_copy_c(umac(n,3), 1,chkdata_edgez(n),1, 1) + call multifab_copy_c(umac_sum(n,3), 1,chkdata_edgez(n),2, 1) + call multifab_copy_c(grad_Epot(n,3), 1,chkdata_edgez(n),3, 1) + end if endif - end if + else ! use_charged_fluid == F + if (advection_type.eq.0) then + call multifab_copy_c(umac(n,1), 1,chkdata_edgex(n),1, 1) + call multifab_copy_c(umac(n,2), 1,chkdata_edgey(n),1, 1) + call multifab_copy_c(umac_sum(n,1), 1,chkdata_edgex(n),2, 1) + call multifab_copy_c(umac_sum(n,2), 1,chkdata_edgey(n),2, 1) + + call multifab_copy_c(mass_fluxes(n,1), 1,chkdata_edgex(n),3, nspecies+1) + call multifab_copy_c(mass_fluxes(n,2), 1,chkdata_edgey(n),3, nspecies+1) + call multifab_copy_c(mass_fluxes_sum(n,1),1,chkdata_edgex(n),(nspecies+1)+3, nspecies+1) + call multifab_copy_c(mass_fluxes_sum(n,2),1,chkdata_edgey(n),(nspecies+1)+3, nspecies+1) + if (dm .eq. 3) then + call multifab_copy_c(umac(n,3), 1,chkdata_edgez(n),1, 1) + call multifab_copy_c(umac_sum(n,3), 1,chkdata_edgez(n),2, 1) + call multifab_copy_c(mass_fluxes(n,3), 1,chkdata_edgez(n),3, nspecies+1) + call multifab_copy_c(mass_fluxes_sum(n,3),1,chkdata_edgez(n),(nspecies+1)+3, nspecies+1) + end if + else ! advection_type \ne 0 + call multifab_copy_c(umac(n,1), 1,chkdata_edgex(n),1, 1) + call multifab_copy_c(umac(n,2), 1,chkdata_edgey(n),1, 1) + call multifab_copy_c(umac_sum(n,1), 1,chkdata_edgex(n),2, 1) + call multifab_copy_c(umac_sum(n,2), 1,chkdata_edgey(n),2, 1) + if (dm .eq. 3) then + call multifab_copy_c(umac(n,3), 1,chkdata_edgez(n),1, 1) + call multifab_copy_c(umac_sum(n,3), 1,chkdata_edgez(n),2, 1) + end if + endif + endif end do do n=1,nlevs diff --git a/staggered_grid/src_lowMach/write_plotfile.f90 b/staggered_grid/src_lowMach/write_plotfile.f90 index 26608473..7fe72c9b 100644 --- a/staggered_grid/src_lowMach/write_plotfile.f90 +++ b/staggered_grid/src_lowMach/write_plotfile.f90 @@ -9,17 +9,18 @@ module write_plotfile_module use eos_check_module use probin_multispecies_module, only: plot_stag, is_nonisothermal use probin_common_module, only: prob_lo, prob_hi, nspecies, plot_base_name, & - algorithm_type, rho0, plot_umac_tavg, plot_Epot_tavg, & - plot_rho_tavg, plot_avg_gradPhiApprox, plot_shifted_vel, & - plot_gradEpot, plot_averaged_vel, plot_debug + algorithm_type, advection_type, rho0, plot_umac_tavg, plot_Epot_tavg, & + plot_rho_tavg, plot_avg_gradPhiApprox, plot_shifted_vel, plot_mass_fluxes, & + plot_mass_fluxes_tavg, plot_gradEpot, plot_averaged_vel, plot_debug, & + plot_charge_fluxes, plot_charge_fluxes_tavg use probin_charged_module, only: use_charged_fluid, electroneutral implicit none contains - subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Epot_avg, & - grad_Epot,gradPhiApprox,istep,dx,time) + subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,mass_fluxes,mass_fluxes_avg,Epot,Epot_avg, & + grad_Epot,gradPhiApprox,charge_fluxes,charge_fluxes_avg,istep,dx,time) type(ml_layout), intent(in) :: mla type(multifab), intent(inout) :: rho(:) @@ -29,10 +30,14 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep type(multifab), intent(in) :: umac(:,:) type(multifab), intent(in) :: umac_avg(:,:) type(multifab), intent(in) :: pres(:) + type(multifab), intent(in) :: mass_fluxes(:,:) + type(multifab), intent(in) :: mass_fluxes_avg(:,:) type(multifab), intent(in) :: Epot(:) type(multifab), intent(in) :: Epot_avg(:) type(multifab), intent(in) :: grad_Epot(:,:) type(multifab), intent(in) :: gradPhiApprox(:,:) + type(multifab), intent(in) :: charge_fluxes(:,:) + type(multifab), intent(in) :: charge_fluxes_avg(:,:) integer, intent(in) :: istep real(kind=dp_t), intent(in) :: dx(:,:),time @@ -73,6 +78,8 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep ! umac_avg averaged :dm (optional) ! umac_avg shifted :dm (optional) ! pressure :1 (pressure) + ! mass_fluxes :(nspecies+1)*dm (optional, default off) + ! mass_fluxes_avg :(nspecies+1)*dm (optional, time-averaged mass fluxes, default off) ! rho_i and pressure nvarsCC = nspecies + 1 @@ -89,13 +96,21 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep if (is_nonisothermal) then nvarsCC = nvarsCC + 1 end if + ! cc velocity if (plot_averaged_vel) then nvarsCC = nvarsCC + dm end if + ! shifted velocity if (plot_shifted_vel) then nvarsCC = nvarsCC + dm end if - ! time-averaged rho + ! mass fluxes--for each individual species this is a vector of dimension dm, plus one vector + ! for the total flux (each individual contribution summed up). This quantity only makes sense + ! advection_type = 0 + if (plot_mass_fluxes.and.(advection_type.eq.0)) then + nvarsCC = nvarsCC + (nspecies+1)*dm + endif + ! time-averaged rho if (plot_rho_tavg) then nvarsCC = nvarsCC + nspecies end if @@ -108,6 +123,12 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep nvarsCC = nvarsCC + dm end if end if + ! time-averaged mass fluxes--for each individual species this is a vector of dimension dm, plus one vector + ! for the total flux (each individual contribution summed up). This quantity only makes sense + ! advection_type = 0 + if (plot_mass_fluxes_tavg.and.(advection_type.eq.0)) then + nvarsCC = nvarsCC + (nspecies+1)*dm + endif if (use_charged_fluid) then ! charge :1 (don't write for electroneutral unless plot_debug=T) @@ -115,6 +136,8 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep ! Epot_avg :1 (optional) ! grad_Epot averaged :dm (optional) ! gradPhiApprox averaged :dm (optional) + ! charge_fluxes :(nspecies+1)*dm + ! tavg_charge_fluxes :(nspecies+1)*dm ! charge if (.not. (electroneutral .and. (.not. plot_debug))) then @@ -134,6 +157,17 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep if (plot_avg_gradPhiApprox) then nvarsCC = nvarsCC + dm end if + ! charge fluxes--for each individual species this is a vector of dimension dm, plus one vector + ! for the total flux (each individual contribution summed up). This quantity only makes sense + ! advection_type = 0 + if (plot_charge_fluxes.and.(advection_type.eq.0)) then + nvarsCC = nvarsCC + (nspecies+1)*dm + endif + ! time-averaged charge fluxes + ! only makes sense for advection_type = 0 + if (plot_charge_fluxes_tavg.and.(advection_type.eq.0)) then + nvarsCC = nvarsCC + (nspecies+1)*dm + endif end if if (boussinesq) then ! Boussinesq @@ -220,6 +254,56 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep plot_names(counter) = "pres" counter = counter + 1 + ! mass fluxes + if (plot_mass_fluxes.and.(advection_type.eq.0)) then + ! for each species + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "mass_flx_x_", n + counter = counter + 1 + enddo + plot_names(counter) = "tot_mass_flx_x" + counter = counter + 1 + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "mass_flx_y_", n + counter = counter + 1 + enddo + plot_names(counter) = "tot_mass_flx_y" + counter = counter + 1 + if (dm > 2) then + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "mass_flx_z_", n + counter = counter + 1 + enddo + plot_names(counter) = "tot_mass_flx_z" + counter = counter + 1 + endif + endif + + ! time-averaged mass fluxes + if (plot_mass_fluxes_tavg.and.(advection_type.eq.0)) then + ! for each species + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "tavg_mass_flx_x_", n + counter = counter + 1 + enddo + plot_names(counter) = "tavg_tot_mass_flx_x" + counter = counter + 1 + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "tavg_mass_flx_y_", n + counter = counter + 1 + enddo + plot_names(counter) = "tavg_tot_mass_flx_y" + counter = counter + 1 + if (dm > 2) then + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "tavg_mass_flx_z_", n + counter = counter + 1 + enddo + plot_names(counter) = "tavg_tot_mass_flx_z" + counter = counter + 1 + endif + endif + if (use_charged_fluid) then if (.not. (electroneutral .and. (.not. plot_debug))) then plot_names(counter) = "charge_density" @@ -256,6 +340,54 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep end if end if + if (plot_charge_fluxes.and.(advection_type.eq.0)) then + ! for each species + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "chrg_flx_x_", n + counter = counter + 1 + enddo + plot_names(counter) = "tot_chrg_flx_x" + counter = counter + 1 + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "chrg_flx_y_", n + counter = counter + 1 + enddo + plot_names(counter) = "tot_chrg_flx_y" + counter = counter + 1 + if (dm > 2) then + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "chrg_flx_z_", n + counter = counter + 1 + enddo + plot_names(counter) = "tot_chrg_flx_z" + counter = counter + 1 + endif + endif + + if (plot_charge_fluxes_tavg.and.(advection_type.eq.0)) then + ! for each species + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "tavg_chrg_flx_x_", n + counter = counter + 1 + enddo + plot_names(counter) = "tavg_tot_chrg_flx_x" + counter = counter + 1 + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "tavg_chrg_flx_y_", n + counter = counter + 1 + enddo + plot_names(counter) = "tavg_tot_chrg_flx_y" + counter = counter + 1 + if (dm > 2) then + do n=1,nspecies + write(plot_names(counter),'(a,i0)') "tavg_chrg_flx_z_", n + counter = counter + 1 + enddo + plot_names(counter) = "tavg_tot_chrg_flx_z" + counter = counter + 1 + endif + endif + end if if (boussinesq) then @@ -395,6 +527,22 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep enddo counter = counter + 1 + ! mass fluxes + if (plot_mass_fluxes.and.(advection_type.eq.0)) then + do i=1,dm + call average_face_to_cc(mla,mass_fluxes(:,i),1,plotdata,counter,nspecies+1) + counter = counter + (nspecies+1) + end do + endif + + ! time-averaged mass fluxes + if (plot_mass_fluxes_tavg.and.(advection_type.eq.0)) then + do i=1,dm + call average_face_to_cc(mla,mass_fluxes_avg(:,i),1,plotdata,counter,nspecies+1) + counter = counter + (nspecies+1) + end do + endif + if (use_charged_fluid) then ! compute total charge, then copy into the correct component @@ -438,6 +586,21 @@ subroutine write_plotfile(mla,rho,rho_avg,rhotot,Temp,umac,umac_avg,pres,Epot,Ep counter = counter + 1 end do end if + + ! charge fluxes + if (plot_charge_fluxes.and.(advection_type.eq.0)) then + do i=1,dm + call average_face_to_cc(mla,charge_fluxes(:,i),1,plotdata,counter,nspecies+1) + counter = counter + (nspecies+1) + end do + endif + ! time-averaged charge fluxes + if (plot_charge_fluxes_tavg.and.(advection_type.eq.0)) then + do i=1,dm + call average_face_to_cc(mla,charge_fluxes_avg(:,i),1,plotdata,counter,nspecies+1) + counter = counter + (nspecies+1) + end do + endif end if if (boussinesq) then diff --git a/staggered_grid/src_multiSpec/LAPACK95/Makefile.Donev b/staggered_grid/src_multiSpec/LAPACK95/Makefile.Donev index 30de0c3f..69ee315d 100644 --- a/staggered_grid/src_multiSpec/LAPACK95/Makefile.Donev +++ b/staggered_grid/src_multiSpec/LAPACK95/Makefile.Donev @@ -1,22 +1,22 @@ -SRC_LAPACK=${HPC_SRC}/FluctHydro/FluctHydro/staggered_grid/src_multiSpec/LAPACK95 +#SRC_LAPACK=~/src/FluctHydro/LowMachFHD/staggered_grid/src_multiSpec/LAPACK95 +SRC_LAPACK=. vpath %.f90 $(SRC_LAPACK) FC = gfortran -ffree-form FC1 = gfortran -ffixed-form -OPTS0 = -O3 +OPTS0 = -O3 -fPIC #MODLIB = -I./lapack95_modules OPTS1 = -c $(OPTS0) OPTS3 = $(OPTS1) $(MODLIB) OPTL = -o OPTLIB = -LAPACK_PATH = /usr/lib64 +LAPACK_PATH = /usr/lib/x86_64-linux-gnu LAPACK95 = lapack95.a -LAPACK77 = $(LAPACK_PATH)/lapack.a -TMG77 = $(LAPACK_PATH)/tmglib.a -BLAS = $(LAPACK_PATH)/blas.a +LAPACK77 = $(LAPACK_PATH)/liblapack.a +BLAS = $(LAPACK_PATH)/libblas64.a include $(SRC_LAPACK)/MakefileGeneric diff --git a/staggered_grid/src_multiSpec/LAPACK95/MakefileGeneric b/staggered_grid/src_multiSpec/LAPACK95/MakefileGeneric index 83e26866..f269e397 100644 --- a/staggered_grid/src_multiSpec/LAPACK95/MakefileGeneric +++ b/staggered_grid/src_multiSpec/LAPACK95/MakefileGeneric @@ -5,6 +5,7 @@ # LIBS = $(LAPACK95) $(TMG77) $(LAPACK77) $(BLAS) SUF = f90 +LIBNAME = liblapack95 XX = 'rm' -f $@; \ 'rm' -f $@.res; \ @@ -17,6 +18,7 @@ YY = $(FC) $(OPTS0) -o $@ $(MODLIB) $@.$(SUF) $(OPTLIB) $(LIBS) .SUFFIXES: .f90 .f .o .$(SUF).o: + $(FC) -v 2>&1 | grep "gcc version" $(FC) $(OPTS3) $< @@ -154,9 +156,9 @@ ZOBJS = la_zgesv.o la_zgesv1.o la_zgesvx.o la_zgesvx1.o \ default: la_auxmod.o single_double single: $(SOBJSS) - rm -f lapack95.a - ar cr lapack95.a $(SOBJSS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(SOBJSS) + ranlib $(LIBNAME).a rm -fr lapack95_modules mkdir lapack95_modules cp *.mod lapack95_modules/ @@ -164,9 +166,9 @@ single: $(SOBJSS) rm -f *_lapack_single.o double: $(DOBJSS) - rm -f lapack95.a - ar cr lapack95.a $(DOBJSS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(DOBJSS) + ranlib $(LIBNAME).a rm -fr lapack95_modules mkdir lapack95_modules cp *.mod lapack95_modules/ @@ -174,9 +176,11 @@ double: $(DOBJSS) rm -f *_lapack_double.o single_double: $(SDOBJS) - rm -f lapack95.a - ar cr lapack95.a $(SDOBJS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(SDOBJS) + ranlib $(LIBNAME).a + echo gcc -shared -fPIC -o $(LIBNAME).so $(SDOBJS) + gcc -shared -fPIC -o $(LIBNAME).so $(SDOBJS) rm -fr lapack95_modules mkdir lapack95_modules cp *.mod lapack95_modules/ @@ -184,9 +188,9 @@ single_double: $(SDOBJS) rm -f *_lapack_single_double.o single_complex: $(SCOBJS) - rm -f lapack95.a - ar cr lapack95.a $(SCOBJS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(SCOBJS) + ranlib $(LIBNAME).a rm -fr lapack95_modules mkdir lapack95_modules cp *.mod lapack95_modules/ @@ -194,9 +198,9 @@ single_complex: $(SCOBJS) rm -f *_lapack_single_complex.o double_dcomplex: $(DZOBJS) - rm -f lapack95.a - ar cr lapack95.a $(DZOBJS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(DZOBJS) + ranlib $(LIBNAME).a rm -fr lapack95_modules mkdir lapack95_modules cp *.mod lapack95_modules/ @@ -204,9 +208,9 @@ double_dcomplex: $(DZOBJS) rm -f *_lapack_double_dcomplex.o single_double_complex: $(SCDOBJS) - rm -f lapack95.a - ar cr lapack95.a $(SCDOBJS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(SCDOBJS) + ranlib $(LIBNAME).a rm -fr lapack95_modules mkdir lapack95_modules cp *.mod lapack95_modules/ @@ -214,9 +218,9 @@ single_double_complex: $(SCDOBJS) rm -f *_lapack_single_double_complex.o single_double_complex_dcomplex: $(SCDZOBJS) - rm -f lapack95.a - ar cr lapack95.a $(SCDZOBJS) - ranlib lapack95.a + rm -f $(LIBNAME).a + ar cr $(LIBNAME).a $(SCDZOBJS) + ranlib $(LIBNAME).a rm -fr lapack95_modules mkdir lapack95_modules 'cp' *.mod lapack95_modules/ diff --git a/staggered_grid/src_multiSpec/compute_mass_fluxdiv.f90 b/staggered_grid/src_multiSpec/compute_mass_fluxdiv.f90 index d17c214b..35fff831 100644 --- a/staggered_grid/src_multiSpec/compute_mass_fluxdiv.f90 +++ b/staggered_grid/src_multiSpec/compute_mass_fluxdiv.f90 @@ -144,7 +144,7 @@ subroutine compute_mass_fluxdiv(mla,rho,rhotot,gradp_baro,Temp, & call electrodiffusive_mass_fluxdiv(mla,rho,Temp,rhoWchi, & diff_mass_flux,diff_mass_fluxdiv, & stoch_mass_flux, & - dx,the_bc_tower, & + dx,stage_time,the_bc_tower, & charge,grad_Epot,Epot, & permittivity,dt,zero_initial_Epot) end if diff --git a/staggered_grid/src_multiSpec/mass_flux_utilities.f90 b/staggered_grid/src_multiSpec/mass_flux_utilities.f90 index 5ef4ba92..eaaa208e 100644 --- a/staggered_grid/src_multiSpec/mass_flux_utilities.f90 +++ b/staggered_grid/src_multiSpec/mass_flux_utilities.f90 @@ -9,7 +9,8 @@ module mass_flux_utilities_module chi_iterations, avg_type, use_multiphase, alpha_gex, n_gex use matrix_utilities use compute_mixture_properties_module - use F95_LAPACK ! Donev: Disabled LAPACK so this builds more easily on different systems + !use F95_LAPACK ! Donev: Disabled LAPACK so this builds more easily on different systems + ! The routines needed from LAPACK are now inside matrix_utilities implicit none @@ -872,12 +873,13 @@ subroutine compute_chi_lapack(nspecies_in,Lambda,chi,W) !========================================================== ! compute chilocal inverse - call dgetrf(nspecies_in, nspecies_in, chilocal, nspecies_in, ipiv, info) - call dgetri(nspecies_in, chilocal, nspecies_in, ipiv, work, nspecies_in, info) - !stop "LAPACK95 dget? disabled" + ! This code uses LAPACK95 but this causes problems so just use old school LAPACK77 instead + !call dgetrf(nspecies_in, nspecies_in, chilocal, nspecies_in, ipiv, info) + !call dgetri(nspecies_in, chilocal, nspecies_in, ipiv, work, nspecies_in, info) + !chi = chilocal ! populate chi with B^(-1) - chi = chilocal + chi = mat_inv(chilocal) case(2) !========================================================== @@ -887,7 +889,8 @@ subroutine compute_chi_lapack(nspecies_in,Lambda,chi,W) ! SVD decomposition of chilocal = U * S * VTranspose; note that chilocal ! is changed. also V=(VT)T, UT = (U)T are needed for pseudoinverse of chilocal. !stop "LAPACK95 la_gesvd disabled" - call la_gesvd(chilocal, S, U, VT) + call mat_svd(chilocal, S, U, VT) + V = transpose(VT) UT = transpose(U) @@ -1507,7 +1510,7 @@ subroutine chol_lapack(sqrtL,nspecies_in) integer :: row, column - call dpotrf_f95(sqrtL,'L', rcond, 'I', info) + call mat_chol(sqrtL,'L', rcond, 'I', info) !stop "LAPACK95 dpotrf_f95 disabled" ! remove all upper-triangular entries and NXN entry that lapack doesn't set to zero diff --git a/staggered_grid/src_multiSpec/matrix_utilities.f90 b/staggered_grid/src_multiSpec/matrix_utilities.f90 index d1b299fd..1b6d1d41 100644 --- a/staggered_grid/src_multiSpec/matrix_utilities.f90 +++ b/staggered_grid/src_multiSpec/matrix_utilities.f90 @@ -6,8 +6,63 @@ module matrix_utilities implicit none private + + integer, parameter, private :: dp=dp_t, wp=dp_t ! Working precision - public :: Dbar2chi_iterative, choldc + public :: Dbar2chi_iterative, choldc, mat_inv, mat_svd, mat_chol + + interface mat_svd + module procedure DGESVD_F95 + end interface + + interface mat_chol + module procedure DPOTRF_F95 + end interface + + INTERFACE ! F77 LAPACK routines needed + + SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, & + & LDVT, WORK, LWORK, INFO ) + IMPORT + CHARACTER(LEN=1), INTENT(IN) :: JOBU, JOBVT + INTEGER, INTENT(IN) :: M, N, LDA, LDU, LDVT, LWORK + INTEGER, INTENT(OUT) :: INFO + REAL(WP), INTENT(OUT) :: S(*) + REAL(WP), INTENT(INOUT) :: A(LDA,*) + REAL(WP), INTENT(OUT) :: U(LDU,*), VT(LDVT,*), WORK(*) + END SUBROUTINE DGESVD + + SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO ) + IMPORT + CHARACTER(LEN=1), INTENT(IN) :: UPLO + INTEGER, INTENT(IN) :: LDA, N + INTEGER, INTENT(OUT) :: INFO + REAL(WP), INTENT(INOUT) :: A(LDA,*) + END SUBROUTINE DPOTRF + + FUNCTION DLANSY( NORM, UPLO, N, A, LDA, WORK ) + IMPORT + REAL(WP) :: DLANSY + CHARACTER(LEN=1), INTENT(IN) :: NORM, UPLO + INTEGER, INTENT(IN) :: LDA, N + REAL(WP), INTENT(IN) :: A( LDA, * ) + REAL(WP), INTENT(OUT) :: WORK( * ) + END FUNCTION DLANSY + + SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, & + & INFO ) + IMPORT + CHARACTER(LEN=1), INTENT(IN) :: UPLO + INTEGER, INTENT(IN) :: LDA, N + INTEGER, INTENT(OUT) :: INFO + REAL(WP), INTENT(IN) :: ANORM + REAL(WP), INTENT(OUT) :: RCOND + INTEGER, INTENT(OUT) :: IWORK( * ) + REAL(WP), INTENT(IN) :: A( LDA, * ) + REAL(WP), INTENT(OUT) :: WORK( * ) + END SUBROUTINE DPOCON + + END INTERFACE contains @@ -19,19 +74,19 @@ module matrix_utilities subroutine Dbar2chi_iterative(nspecies_local,num_iterations,D_bar,Xk,molmass_local,chi) integer, intent(in) :: nspecies_local integer, intent(in) :: num_iterations - real(kind=dp_t), intent(in) :: D_bar(1:nspecies_local,1:nspecies_local) - real(kind=dp_t), intent(in) :: Xk(1:nspecies_local), molmass_local(1:nspecies_local) - real(kind=dp_t), intent(out) :: chi(1:nspecies_local,1:nspecies_local) + real(kind=wp), intent(in) :: D_bar(1:nspecies_local,1:nspecies_local) + real(kind=wp), intent(in) :: Xk(1:nspecies_local), molmass_local(1:nspecies_local) + real(kind=wp), intent(out) :: chi(1:nspecies_local,1:nspecies_local) ! Local variables - real(kind=dp_t) :: term1, term2, MWmix - real(kind=dp_t) :: Di(1:nspecies_local) - real(kind=dp_t) :: Deltamat(1:nspecies_local,1:nspecies_local), Zmat(1:nspecies_local,1:nspecies_local) - real(kind=dp_t), dimension(1:nspecies_local,1:nspecies_local) :: Pmat, Jmat - real(kind=dp_t), dimension(1:nspecies_local) :: Minv, Mmat - real(kind=dp_t), dimension(1:nspecies_local,1:nspecies_local) :: PJ, matrix1, matrix2 - real(kind=dp_t) :: scr - real(kind=dp_t) :: Ykp(1:nspecies_local), Xkp(1:nspecies_local) + real(kind=wp) :: term1, term2, MWmix + real(kind=wp) :: Di(1:nspecies_local) + real(kind=wp) :: Deltamat(1:nspecies_local,1:nspecies_local), Zmat(1:nspecies_local,1:nspecies_local) + real(kind=wp), dimension(1:nspecies_local,1:nspecies_local) :: Pmat, Jmat + real(kind=wp), dimension(1:nspecies_local) :: Minv, Mmat + real(kind=wp), dimension(1:nspecies_local,1:nspecies_local) :: PJ, matrix1, matrix2 + real(kind=wp) :: scr + real(kind=wp) :: Ykp(1:nspecies_local), Xkp(1:nspecies_local) integer :: i, j, k, ii, jj @@ -161,15 +216,15 @@ subroutine Dbar2chi_iterative(nspecies_local,num_iterations,D_bar,Xk,molmass_loc ! upon return the lower triangle and diagonal are overwritten by the cholesky factor subroutine choldc(a,np) integer :: np - real(kind=dp_t), intent(inout) :: a(np,np) + real(kind=wp), intent(inout) :: a(np,np) - real(kind=dp_t) :: p(np), dij(np,np) - real(kind=dp_t) :: dd(np,np) - real(kind=dp_t) :: yy(np), mwmix + real(kind=wp) :: p(np), dij(np,np) + real(kind=wp) :: dd(np,np) + real(kind=wp) :: yy(np), mwmix integer :: i, j, k, ii, jj - real(kind=dp_t) :: sum1 - real(kind=dp_t) :: small_number = 0.0d0 ! Some tolerance + real(kind=wp) :: sum1 + real(kind=wp) :: small_number = 0.0d0 ! Some tolerance integer :: idiag,ising @@ -240,4 +295,610 @@ subroutine choldc(a,np) end subroutine +!------------------------------------------------- +! Code adopted by Aleks Donev from +! https://fortranwiki.org/fortran/show/Matrix+inversion +!------------------------------------------------- + +function mat_inv(A) result(Ainv) + real(wp), dimension(:,:), intent(in) :: A + real(wp), dimension(size(A,1),size(A,2)) :: Ainv + + integer :: n + + if(size(A,1)/=size(A,2)) then + stop "Matrix must be square" + end if + + n = size(A,1) + select case(n) + case(1) + Ainv(1,1)=1.0_wp/A(1,1) + case(2) + Ainv=matinv2(A) + case(3) + Ainv=matinv3(A) + case(4) + Ainv=matinv4(A) + case default + Ainv=matinv2(A) + end select + +end function + +! Returns the inverse of a matrix calculated by finding the LU +! decomposition. Depends on LAPACK. +function inv(A) result(Ainv_wp) + real(wp), dimension(:,:), intent(in) :: A + real(wp), dimension(size(A,1),size(A,2)) :: Ainv_wp + + real(dp), dimension(size(A,1),size(A,2)) :: Ainv ! Use LAPACK double precision routines + + real(dp), dimension(size(A,1)) :: work ! work array for LAPACK + integer, dimension(size(A,1)) :: ipiv ! pivot indices + integer :: n, info + + ! External procedures defined in LAPACK + external DGETRF + external DGETRI + + ! Store A in Ainv to prevent it from being overwritten by LAPACK + Ainv = A + n = size(A,1) + + ! DGETRF computes an LU factorization of a general M-by-N matrix A + ! using partial pivoting with row interchanges. + call DGETRF(n, n, Ainv, n, ipiv, info) + + if (info /= 0) then + stop 'Matrix is numerically singular!' + end if + + ! DGETRI computes the inverse of a matrix using the LU factorization + ! computed by DGETRF. + call DGETRI(n, Ainv, n, ipiv, work, n, info) + + if (info /= 0) then + stop 'Matrix inversion failed!' + end if + + Ainv_wp = Ainv ! Convert precisions if needed + +end function inv + + pure function matinv2(A) result(B) + !! Performs a direct calculation of the inverse of a 2×2 matrix. + real(wp), intent(in) :: A(2,2) !! Matrix + real(wp) :: B(2,2) !! Inverse matrix + real(wp) :: detinv + + ! Calculate the inverse determinant of the matrix + detinv = 1/(A(1,1)*A(2,2) - A(1,2)*A(2,1)) + + ! Calculate the inverse of the matrix + B(1,1) = +detinv * A(2,2) + B(2,1) = -detinv * A(2,1) + B(1,2) = -detinv * A(1,2) + B(2,2) = +detinv * A(1,1) + end function + + pure function matinv3(A) result(B) + !! Performs a direct calculation of the inverse of a 3×3 matrix. + real(wp), intent(in) :: A(3,3) !! Matrix + real(wp) :: B(3,3) !! Inverse matrix + real(wp) :: detinv + + ! Calculate the inverse determinant of the matrix + detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + ! Calculate the inverse of the matrix + B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2)) + B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1)) + B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1)) + B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2)) + B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1)) + B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1)) + B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2)) + B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1)) + B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1)) + end function + + pure function matinv4(A) result(B) + !! Performs a direct calculation of the inverse of a 4×4 matrix. + real(wp), intent(in) :: A(4,4) !! Matrix + real(wp) :: B(4,4) !! Inverse matrix + real(wp) :: detinv + + ! Calculate the inverse determinant of the matrix + detinv = & + 1/(A(1,1)*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2)))& + - A(1,2)*(A(2,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1)))& + + A(1,3)*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))& + - A(1,4)*(A(2,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(2,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(2,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))) + + ! Calculate the inverse of the matrix + B(1,1) = detinv*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))) + B(2,1) = detinv*(A(2,1)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(2,3)*(A(3,1)*A(4,4)-A(3,4)*A(4,1))+A(2,4)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))) + B(3,1) = detinv*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))) + B(4,1) = detinv*(A(2,1)*(A(3,3)*A(4,2)-A(3,2)*A(4,3))+A(2,2)*(A(3,1)*A(4,3)-A(3,3)*A(4,1))+A(2,3)*(A(3,2)*A(4,1)-A(3,1)*A(4,2))) + B(1,2) = detinv*(A(1,2)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(1,3)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(1,4)*(A(3,3)*A(4,2)-A(3,2)*A(4,3))) + B(2,2) = detinv*(A(1,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(1,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(1,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1))) + B(3,2) = detinv*(A(1,1)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(1,2)*(A(3,1)*A(4,4)-A(3,4)*A(4,1))+A(1,4)*(A(3,2)*A(4,1)-A(3,1)*A(4,2))) + B(4,2) = detinv*(A(1,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(1,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(1,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))) + B(1,3) = detinv*(A(1,2)*(A(2,3)*A(4,4)-A(2,4)*A(4,3))+A(1,3)*(A(2,4)*A(4,2)-A(2,2)*A(4,4))+A(1,4)*(A(2,2)*A(4,3)-A(2,3)*A(4,2))) + B(2,3) = detinv*(A(1,1)*(A(2,4)*A(4,3)-A(2,3)*A(4,4))+A(1,3)*(A(2,1)*A(4,4)-A(2,4)*A(4,1))+A(1,4)*(A(2,3)*A(4,1)-A(2,1)*A(4,3))) + B(3,3) = detinv*(A(1,1)*(A(2,2)*A(4,4)-A(2,4)*A(4,2))+A(1,2)*(A(2,4)*A(4,1)-A(2,1)*A(4,4))+A(1,4)*(A(2,1)*A(4,2)-A(2,2)*A(4,1))) + B(4,3) = detinv*(A(1,1)*(A(2,3)*A(4,2)-A(2,2)*A(4,3))+A(1,2)*(A(2,1)*A(4,3)-A(2,3)*A(4,1))+A(1,3)*(A(2,2)*A(4,1)-A(2,1)*A(4,2))) + B(1,4) = detinv*(A(1,2)*(A(2,4)*A(3,3)-A(2,3)*A(3,4))+A(1,3)*(A(2,2)*A(3,4)-A(2,4)*A(3,2))+A(1,4)*(A(2,3)*A(3,2)-A(2,2)*A(3,3))) + B(2,4) = detinv*(A(1,1)*(A(2,3)*A(3,4)-A(2,4)*A(3,3))+A(1,3)*(A(2,4)*A(3,1)-A(2,1)*A(3,4))+A(1,4)*(A(2,1)*A(3,3)-A(2,3)*A(3,1))) + B(3,4) = detinv*(A(1,1)*(A(2,4)*A(3,2)-A(2,2)*A(3,4))+A(1,2)*(A(2,1)*A(3,4)-A(2,4)*A(3,1))+A(1,4)*(A(2,2)*A(3,1)-A(2,1)*A(3,2))) + B(4,4) = detinv*(A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))) + end function + +!----------------------------------------------- +! A. Donev adopted from LAPACK90 interface +!----------------------------------------------- +SUBROUTINE DGESVD_F95( A, S, U, VT, WW, JOB, INFO ) +! +! -- LAPACK95 interface driver routine (version 3.0) -- +! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK +! September, 2000 +! +! .. USE STATEMENTS .. + !USE LA_PRECISION, ONLY: WP => DP + !USE LA_AUXMOD, ONLY: ERINFO, LSAME +! .. IMPLICIT STATEMENT .. + IMPLICIT NONE +! .. SCALAR ARGUMENTS .. + CHARACTER(LEN=1), OPTIONAL, INTENT(IN) :: JOB + INTEGER, INTENT(OUT), OPTIONAL :: INFO +! .. ARRAY ARGUMENTS .. + REAL(WP), INTENT(INOUT) :: A(:,:) + REAL(WP), INTENT(OUT) :: S(:) + REAL(WP), INTENT(OUT), OPTIONAL :: WW(:) + REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: U(:,:), VT(:,:) +!---------------------------------------------------------------------- +! +! Purpose +! ======= +! +! LA_GESVD and LA_GESDD compute the singular values and, +! optionally, the left and/or right singular vectors from the singular +! value decomposition (SVD) of a real or complex m by n matrix A. The +! SVD of A is written +! A = U * SIGMA * V^H +! where SIGMA is an m by n matrix which is zero except for its +! min(m, n) diagonal elements, U is an m by m orthogonal (unitary) +! matrix, and V is an n by n orthogonal (unitary) matrix. The diagonal +! elements of SIGMA , i.e., the values +! +! sigma(i)= SIGMA(i,i), i = 1, 2,..., min(m, n) +! are the singular values of A; they are real and non-negative, and are +! returned in descending order. The first min(m, n) columns of U and V +! are the left and right singular vectors of A, respectively. +! LA_GESDD solves the same problem as LA_GESVD but uses a divide and +! conquer method if singular vectors are desired. For large matrices it +! is usually much faster than LA_GESVD when singular vectors are +! desired, but uses more workspace. +! +! Note: The routine returns V^H , not V . +! +! ======== +! +! SUBROUTINE LA_GESVD / LA_GESDD( A, S, U=u, VT=vt, & +! WW=ww, JOB=job, INFO=info ) +! (), INTENT(INOUT) :: A(:,:) +! REAL(), INTENT(OUT) :: S(:) +! (), INTENT(OUT), OPTIONAL :: U(:,:), VT(:,:) +! REAL(), INTENT(OUT), OPTIONAL :: WW(:) +! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: JOB +! INTEGER, INTENT(OUT), OPTIONAL :: INFO +! where +! ::= REAL | COMPLEX +! ::= KIND(1.0) | KIND(1.0D0) +! +! Arguments +! ========= +! +! A (input/output) REAL or COMPLEX array, shape (:, :) with +! size(A, 1) = m and size(A, 2) = n. +! On entry, the matrix A. +! On exit, if JOB = 'U' and U is not present, then A is +! overwritten with the first min(m, n) columns of U (the left +! singular vectors, stored columnwise). +! If JOB = 'V' and VT is not present, then A is overwritten with +! the first min(m, n) rows of V^H (the right singular vectors, +! stored rowwise). +! In all cases the original contents of A are destroyed. +! S (output) REAL array, shape (:) with size(S) = min(m, n). +! The singular values of A, sorted so that S(i) >= S(i+1). +! U Optional (output) REAL or COMPLEX array, shape (:, :) with +! size(U, 1) = m and size(U, 2) = m or min(m, n). +! If size(U, 2) = m, U contains the m by m matrix U . +! If size(U; 2) = min(m, n), U contains the first min(m, n) +! columns of U (the left singular vectors, stored columnwise). +! VT Optional (output) REAL or COMPLEX array, shape (:, :) with +! size(VT, 1) = n or min(m, n) and size(VT, 2) = n. +! If size(VT, 1) = n , VT contains the n by n matrix V^H . +! If size(VT, 1) = min(m, n), VT contains the first min(m, n) +! rows of V^H (the right singular vectors, stored rowwise). +! WW Optional (output) REAL array, shape (:) with size(WW) = +! min(m, n) - 1 +! If INFO > 0, WW contains the unconverged superdiagonal elements +! of an upper bidiagonal matrix B whose diagonal is in SIGMA (not +! necessarily sorted). B has the same singular values as A. +! Note: WW is a dummy argument for LA_GESDD. +! JOB Optional (input) CHARACTER(LEN=1). +! = 'N': neither columns of U nor rows of V^H are returned in +! array A. +! = 'U': if U is not present, the first min(m, n) columns of U +! (the left singular vectors) are returned in array A; +! = 'V': if VT is not present, the first min(m, n) rows of V^H +! (the right singular vectors) are returned in array A; +! Default value: 'N'. +! INFO Optional (output) INTEGER. +! = 0: successful exit. +! < 0: if INFO = -i, the i-th argument had an illegal value. +! > 0: The algorithm did not converge. +! If INFO is not present and an error occurs, then the program is +! terminated with an error message. +!---------------------------------------------------------------------- +! .. LOCAL PARAMETERS .. + CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GESVD' +! .. LOCAL SCALARS .. + CHARACTER(LEN=1) :: LJOB + CHARACTER(LEN=1) :: LJOBU, LJOBVT + INTEGER, SAVE :: LWORK = 0 + INTEGER :: N, M, LINFO, LD, ISTAT, ISTAT1, S1U, S2U, S1VT, S2VT, & + NN, MN, SWW +! .. LOCAL ARRAYS .. + REAL(WP), TARGET :: LLU(1,1), LLVT(1,1) + REAL(WP), POINTER :: WORK(:) +! .. INTRINSIC FUNCTIONS .. + INTRINSIC MIN, MAX, PRESENT, SIZE +! .. EXECUTABLE STATEMENTS .. + LINFO = 0; ISTAT = 0; M = SIZE(A,1); N = SIZE(A,2) + LD = MAX(1,M); MN = MIN(M,N) + IF( PRESENT(JOB) )THEN; LJOB = JOB; ELSE; LJOB = 'N'; ENDIF + IF( PRESENT(U) )THEN; S1U = SIZE(U,1); S2U = SIZE(U,2) + ELSE; S1U = 1; S2U = 1; END IF + IF( PRESENT(VT) )THEN; S1VT = SIZE(VT,1); S2VT = SIZE(VT,2) + ELSE; S1VT = 1; S2VT = 1; END IF + IF( PRESENT(WW) )THEN; SWW = SIZE(WW); ELSE; SWW = MN-1; ENDIF +! .. TEST THE ARGUMENTS + IF( M < 0 .OR. N < 0 )THEN; LINFO = -1 + ELSE IF( SIZE( S ) /= MN )THEN; LINFO = -2 + ELSE IF( PRESENT(U) .AND. ( S1U /= M .OR. & + ( S2U /= M .AND. S2U /= MN ) ) )THEN; LINFO = -3 + ELSE IF( PRESENT(VT) .AND. ( ( S1VT /= N .AND. S1VT /= MN ) & + .OR. S2VT /= N ) )THEN; LINFO = -4 + ELSE IF( SWW /= MN-1 .AND. MN > 0 ) THEN; LINFO = -5 + ELSE IF( PRESENT(JOB) .AND. ( .NOT. ( LSAME(LJOB,'U') .OR. & + LSAME(LJOB,'V') .OR. LSAME(LJOB,'N') ) .OR. & + LSAME(LJOB,'U') .AND. PRESENT(U) .OR. & + LSAME(LJOB,'V') .AND. PRESENT(VT)) )THEN; LINFO = -6 + ELSE + IF( PRESENT(U) )THEN + IF( S2U == M )THEN; LJOBU = 'A'; ELSE; LJOBU = 'S'; ENDIF + ELSE; IF( LSAME(LJOB,'U') ) THEN; LJOBU = 'O' + ELSE; LJOBU = 'N'; ENDIF + ENDIF + IF( PRESENT(VT) )THEN + IF( S1VT == N )THEN; LJOBVT = 'A'; ELSE; LJOBVT = 'S'; ENDIF + ELSE; IF( LSAME(LJOB,'V') )THEN; LJOBVT = 'O' + ELSE; LJOBVT = 'N'; ENDIF + ENDIF + IF( ISTAT == 0 )THEN + NN = MAX( 5, 3*MN + MAX(M,N), 5*MN ) + LWORK = MAX( 1, NN, LWORK ); ALLOCATE(WORK(LWORK), STAT=ISTAT) + IF( ISTAT /= 0 )THEN; DEALLOCATE(WORK,STAT=ISTAT1) + LWORK = MAX( 1, NN); ALLOCATE(WORK(LWORK), STAT=ISTAT) + IF( ISTAT == 0) CALL ERINFO( -200, SRNAME, LINFO ) + END IF + END IF + IF( ISTAT == 0 ) THEN + IF( PRESENT(U) ) THEN + IF ( PRESENT(VT) )THEN + CALL DGESVD( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), & + VT, MAX(1,S1VT), WORK, LWORK, LINFO ) + ELSE + CALL DGESVD( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), & + & LLVT, MAX(1,S1VT), WORK, LWORK, LINFO ) + ENDIF + ELSE + IF ( PRESENT(VT) )THEN + CALL DGESVD( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), & + & VT, MAX(1,S1VT), WORK, LWORK, LINFO ) + ELSE + CALL DGESVD( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), & + & LLVT, MAX(1,S1VT), WORK, LWORK, LINFO ) + ENDIF + ENDIF + LWORK = INT(WORK(1)+1) + IF( LINFO > 0 .AND. PRESENT(WW) ) WW(1:MN-1) = WORK(2:MN) + ELSE; LINFO = -100; ENDIF + DEALLOCATE(WORK, STAT=ISTAT1) + ENDIF + CALL ERINFO(LINFO,SRNAME,INFO,ISTAT) +END SUBROUTINE DGESVD_F95 + +SUBROUTINE DPOTRF_F95( A, UPLO, RCOND, NORM, INFO ) +! +! -- LAPACK95 interface driver routine (version 3.0) -- +! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK +! September, 2000 +! +! .. USE STATEMENTS .. + !USE LA_PRECISION, ONLY: WP => DP + !USE LA_AUXMOD, ONLY: ERINFO, LSAME +! .. IMPLICIT STATEMENT .. + IMPLICIT NONE +! .. CHARACTER ARGUMENTS .. + CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM, UPLO +! .. SCALAR ARGUMENTS .. + INTEGER, INTENT(OUT), OPTIONAL :: INFO + REAL(WP), INTENT(OUT), OPTIONAL :: RCOND +! .. ARRAY ARGUMENTS .. + REAL(WP), INTENT(INOUT) :: A(:,:) +!----------------------------------------------------------------- +! +! Purpose +! ======= +! +! LA_POTRF computes the Cholesky factorization of a real symmetric or +! complex Hermitian positive definite matrix A. +! +! The factorization has the form +! A = U**H * U, if UPLO = 'U', or +! A = L * L**H, if UPLO = 'L', +! where U is an upper triangular matrix and L is lower triangular. +! +! This is the block version of the algorithm, calling Level 3 BLAS. +! +! LA_POTRF optionally estimates the reciprocal of the condition number +! (in the 1-norm) of a real symmetric or complex Hermitian positive +! definite matrix A. +! An estimate is obtained for norm(inv(A)), and the reciprocal of the +! condition number is computed as RCOND = 1 / (norm(A) * norm(inv(A))). +! +! ======= +! +! SUBROUTINE LA_POTRF( A, UPLO, RCOND, NORM, INFO ) +! (), INTENT(INOUT) :: A(:,:) +! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO +! REAL(), INTENT(OUT), OPTIONAL :: RCOND +! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM +! INTEGER, INTENT(OUT), OPTIONAL :: INFO +! where +! ::= REAL | COMPLEX +! ::= KIND(1.0) | KIND(1.0D0) +! +! Defaults +! ======== +! +! 1. If UPLO is not present then UPLO = 'U' is assumed. +! +! Arguments +! ========= +! +! A (input/output) either REAL or COMPLEX square array, +! shape (:,:), size(A,1) == size(A,2) >= 0. +! On entry, the symmetric (Hermitian) matrix A. +! If UPLO = 'U', the upper triangular part of A contains +! the upper triangular part of the matrix A, and the +! strictly lower triangular part of A is not referenced. +! If UPLO = 'L', the lower triangular part of A contains +! the lower triangular part of the matrix A, and the +! strictly upper triangular part of A is not referenced. +! On exit, if INFO = 0, the factor U or L from the Cholesky +! factorization A = U**H*U or A = L*L**H. +! +! UPLO Optional, (input) CHARACTER*1 +! If UPLO is present then: +! = 'U': Upper triangle of A is stored; +! = 'L': Lower triangle of A is stored. +! otherwise UPLO = 'U' is assumed. +! +! RCOND Optional (output) REAL +! The reciprocal of the condition number of the matrix A +! computed as RCOND = 1/(norm(A) * norm(inv(A))). +! NORM Optional (input) CHARACTER*1 +! Specifies whether the 1-norm condition number or the +! infinity-norm condition number is required: +! If NORM is present then: +! = '1', 'O' or 'o': 1-norm; +! = 'I' or 'i': infinity-norm. +! otherwise NORM = '1' is used. +! +! INFO Optional, (output) INTEGER +! If INFO is present: +! = 0: successful exit +! < 0: if INFO = -i, the i-th argument had an illegal value +! > 0: if INFO = i, the leading minor of order i is not +! positive definite, and the factorization could not be +! completed. +! If INFO is not present and an error occurs, then the program +! is terminated with an error message. +! +! -------------------------------------- +! .. LOCAL PARAMETERS .. + CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_POTRF' +! .. LOCAL SCALARS .. + CHARACTER(LEN=1) :: LNORM, LUPLO + INTEGER :: LINFO, N, ISTAT, ISTAT1, LD + REAL(WP) :: ANORM +! .. LOCAL POINTERS .. + INTEGER, POINTER :: IWORK(:) + REAL(WP), POINTER :: WORK(:) +! .. INTRINSIC FUNCTIONS .. + INTRINSIC PRESENT, MAX +! .. EXECUTABLE STATEMENTS .. + LINFO = 0; N = SIZE(A,1); LD = MAX(1,N); ISTAT = 0 + IF( PRESENT(UPLO) ) THEN; LUPLO = UPLO; ELSE; LUPLO = 'U'; END IF + IF( PRESENT(NORM) ) THEN; LNORM = NORM; ELSE; LNORM = '1'; END IF +! .. TEST THE ARGUMENTS + IF( SIZE( A, 2 ) /= N .AND. N < 0 )THEN; LINFO = -1 + ELSE IF( .NOT.LSAME(LUPLO,'U') .AND. .NOT.LSAME(LUPLO,'L') )THEN; LINFO = -2 + ELSE IF( ( .NOT.PRESENT(RCOND) .AND. PRESENT(NORM) ) .OR. & + ( .NOT.LSAME(LNORM,'I') .AND. .NOT.LSAME(LNORM,'O') & + .AND. LNORM /= '1' ) ) THEN; LINFO = -4 + ELSE IF( N > 0 )THEN + IF( PRESENT(RCOND) ) THEN +! .. COMPUTE THE NORM OF THE MATRIX A + ALLOCATE(WORK(N), STAT=ISTAT) + IF( ISTAT == 0 )THEN; ANORM = DLANSY( LNORM, LUPLO, LD, A, N, WORK ) + ELSE; LINFO = -100; END IF + DEALLOCATE(WORK, STAT=ISTAT1) + END IF +! + IF( LINFO == 0 ) THEN +! .. COMPUTE THE CHOLESKY FACTORS OF THE MATRIX A + CALL DPOTRF( LUPLO, N, A, LD, LINFO ) +! + IF( PRESENT(RCOND) .AND. LINFO == 0 ) THEN +! .. COMPUTE THE RECIPROCAL OF THE CONDITION NUMBER OF A + IF( ANORM == 0.0_WP )THEN; RCOND = 0.0_WP + ELSE; ALLOCATE(WORK(3*N), IWORK(N), STAT=ISTAT) + IF( ISTAT == 0 )THEN + CALL DPOCON( LUPLO, N, A, LD, ANORM, RCOND, & + WORK, IWORK, LINFO ) + ELSE; LINFO = -100; END IF + DEALLOCATE(WORK, IWORK, STAT=ISTAT1) + END IF + END IF + END IF + ELSE IF( PRESENT(RCOND) ) THEN; RCOND = 1.0_WP; ENDIF + CALL ERINFO(LINFO,SRNAME,INFO,ISTAT) +END SUBROUTINE DPOTRF_F95 + + +!------------------------------------------ + LOGICAL FUNCTION LSAME( CA, CB ) +! +! PURPOSE +! ======= +! +! LSAME TESTS IF CA IS THE SAME LETTER AS CB REGARDLESS OF CASE. +! +! PARAMETERS +! ========== +! +! CA (INPUT) CHARACTER*1 +! CB (INPUT) CHARACTER*1 +! CHARACTERS TO BE COMPARED. +! +! .. SCALAR ARGUMENTS .. + CHARACTER*1, INTENT(IN) :: CA, CB +! .. PARAMETERS .. + INTEGER, PARAMETER :: IOFF=32 +! .. LOCAL SCALARS .. + INTEGER :: INTA, INTB, ZCODE +! .. INTRINSIC FUNCTIONS .. + INTRINSIC ICHAR +! +! .. EXECUTABLE STATEMENTS .. +! +! TEST IF THE CHARACTERS ARE EQUAL +! + LSAME = CA == CB +! +! NOW TEST FOR EQUIVALENCE +! + IF( .NOT.LSAME )THEN +! +! USE 'Z' RATHER THAN 'A' SO THAT ASCII CAN BE DETECTED ON PRIME +! MACHINES, ON WHICH ICHAR RETURNS A VALUE WITH BIT 8 SET. +! ICHAR('A') ON PRIME MACHINES RETURNS 193 WHICH IS THE SAME AS +! ICHAR('A') ON AN EBCDIC MACHINE. +! + ZCODE = ICHAR( 'Z' ) +! + INTA = ICHAR( CA ) + INTB = ICHAR( CB ) +! + IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 )THEN +! +! ASCII IS ASSUMED - ZCODE IS THE ASCII CODE OF EITHER LOWER OR +! UPPER CASE 'Z'. +! + IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 + IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 +! + ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 )THEN +! +! EBCDIC IS ASSUMED - ZCODE IS THE EBCDIC CODE OF EITHER LOWER OR +! UPPER CASE 'Z'. +! + IF( INTA.GE.129 .AND. INTA.LE.137 .OR. & +! & INTA.GE.145 .AND. INTA.LE.153 .OR. & + & INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 + IF( INTB.GE.129 .AND. INTB.LE.137 .OR. & + & INTB.GE.145 .AND. INTB.LE.153 .OR. & + & INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 +! + ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 )THEN +! +! ASCII IS ASSUMED, ON PRIME MACHINES - ZCODE IS THE ASCII CODE +! PLUS 128 OF EITHER LOWER OR UPPER CASE 'Z'. +! + IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 + IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 + ENDIF + LSAME = INTA == INTB + ENDIF + END FUNCTION LSAME + +!------------------------------------------------ + SUBROUTINE ERINFO(LINFO, SRNAME, INFO, ISTAT) +! +! -- LAPACK95 interface driver routine (version 3.0) -- +! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK +! September, 2000 +! +! .. IMPLICIT STATEMENT .. + IMPLICIT NONE +! .. SCALAR ARGUMENTS .. + CHARACTER( LEN = * ), INTENT(IN) :: SRNAME + INTEGER , INTENT(IN) :: LINFO + INTEGER , INTENT(OUT), OPTIONAL :: INFO + INTEGER , INTENT(IN), OPTIONAL :: ISTAT +! .. EXECUTABLE STATEMENTS .. +! IF( ( LINFO < 0 .AND. LINFO > -200 ) .OR. & +! & ( LINFO > 0 .AND. .NOT.PRESENT(INFO) ) )THEN + IF( ( ( LINFO < 0 .AND. LINFO > -200 ) .OR. LINFO > 0 ) & + & .AND. .NOT.PRESENT(INFO) )THEN + WRITE (*,*) 'Program terminated in LAPACK95 subroutine ',SRNAME + WRITE (*,*) 'Error indicator, INFO = ',LINFO + IF( PRESENT(ISTAT) )THEN + IF( ISTAT /= 0 ) THEN + IF( LINFO == -100 )THEN + WRITE (*,*) 'The statement ALLOCATE causes STATUS = ', & + & ISTAT + ELSE + WRITE (*,*) 'LINFO = ', LINFO, ' not expected' + END IF + END IF + END IF + STOP + ELSE IF( LINFO <= -200 ) THEN + WRITE(*,*) '++++++++++++++++++++++++++++++++++++++++++++++++' + WRITE(*,*) '*** WARNING, INFO = ', LINFO, ' WARNING ***' + IF( LINFO == -200 )THEN + WRITE(*,*) & + & 'Could not allocate sufficient workspace for the optimum' + WRITE(*,*) & + & 'blocksize, hence the routine may not have performed as' + WRITE(*,*) 'efficiently as possible' + ELSE + WRITE(*,*) 'Unexpected warning' + END IF + WRITE(*,*) '++++++++++++++++++++++++++++++++++++++++++++++++' + END IF + IF( PRESENT(INFO) ) THEN + INFO = LINFO + END IF + END SUBROUTINE ERINFO + end module diff --git a/staggered_grid/src_multiSpec/probin_multispecies.f90 b/staggered_grid/src_multiSpec/probin_multispecies.f90 index 79093f63..885cd575 100644 --- a/staggered_grid/src_multiSpec/probin_multispecies.f90 +++ b/staggered_grid/src_multiSpec/probin_multispecies.f90 @@ -26,6 +26,7 @@ module probin_multispecies_module integer :: midpoint_stoch_mass_flux_type integer :: avg_type integer :: mixture_type + logical :: additive_noise real(kind=dp_t) :: rhotot_bc(3,2) ! used for boussinesq algorithm - not part of namelist ! code sets this in src_lowMach/define_bc_tower.f90 @@ -81,6 +82,8 @@ module probin_multispecies_module ! See compute_mixture_properties.f90 for values supported at present ! The default mixture_type=0 means no dependence on composition + namelist /probin_multispecies/ additive_noise ! additive or multiplicative noise + contains subroutine probin_multispecies_init() @@ -126,6 +129,7 @@ subroutine probin_multispecies_init() midpoint_stoch_mass_flux_type = 1 avg_type = 1 mixture_type = 0 + additive_noise = .false. ! read from input file need_inputs = .true. @@ -268,6 +272,11 @@ subroutine probin_multispecies_init() call get_command_argument(farg, value = fname) read(fname, *) mixture_type + case ('--additive_noise') + farg = farg + 1 + call get_command_argument(farg, value = fname) + read(fname, *) additive_noise + case ('--') farg = farg + 1 exit diff --git a/staggered_grid/src_multiSpec/stochastic_mass_fluxdiv.f90 b/staggered_grid/src_multiSpec/stochastic_mass_fluxdiv.f90 index cbf83528..5907ee76 100644 --- a/staggered_grid/src_multiSpec/stochastic_mass_fluxdiv.f90 +++ b/staggered_grid/src_multiSpec/stochastic_mass_fluxdiv.f90 @@ -1,4 +1,4 @@ -module stochastic_mass_fluxdiv_module + module stochastic_mass_fluxdiv_module use bl_types use bl_space @@ -17,19 +17,23 @@ module stochastic_mass_fluxdiv_module use BoxLibRNGs use bl_rng_module use probin_common_module, only: k_B, molmass, variance_coef_mass, use_bl_rng, nspecies - use probin_multispecies_module, only: correct_flux, is_nonisothermal + use probin_multispecies_module, only: correct_flux, is_nonisothermal, additive_noise implicit none private public :: init_mass_stochastic, destroy_mass_stochastic, & - stochastic_mass_fluxdiv, fill_mass_stochastic + stochastic_mass_fluxdiv, fill_mass_stochastic, & + init_additive_sqrtLonsager, destroy_additive_sqrtLonsager, & + additive_sqrtLonsager_fc ! stochastic fluxes for mass densities are face-centered type(multifab), allocatable, save :: stoch_W_fc(:,:,:) + type(multifab), allocatable, save :: additive_sqrtLonsager_fc(:,:) + integer, save :: n_rngs ! how many random number stages contains @@ -82,7 +86,11 @@ subroutine stochastic_mass_fluxdiv(mla,rho,rhotot,sqrtLonsager_fc, & ! compute variance X sqrtLonsager-face X W(0,1) do n=1,nlevs do i=1,dm - call matvec_mul(mla, stoch_mass_flux(n,i), sqrtLonsager_fc(n,i), nspecies) + if (additive_noise) then + call matvec_mul(mla, stoch_mass_flux(n,i), additive_sqrtLonsager_fc(n,i), nspecies) + else + call matvec_mul(mla, stoch_mass_flux(n,i), sqrtLonsager_fc(n,i), nspecies) + end if call multifab_mult_mult_s(stoch_mass_flux(n,i), variance, 0) end do end do @@ -366,5 +374,45 @@ subroutine stoch_mass_bc_3d(sflux,ng_f,idim,lo,hi,phys_bc) end if end subroutine stoch_mass_bc_3d + + subroutine init_additive_sqrtLonsager(mla) + + type(ml_layout), intent(in ) :: mla + + ! local + integer n,nlevs,i,dm + + nlevs = mla%nlevel + dm = mla%dim + + allocate(additive_sqrtLonsager_fc(nlevs,dm)) + + do n=1,nlevs + do i=1,dm + call multifab_build_edge(additive_sqrtLonsager_fc(n,i), mla%la(n), nspecies**2, 0, i) + end do + end do + + end subroutine init_additive_sqrtLonsager + + subroutine destroy_additive_sqrtLonsager(mla) + + type(ml_layout), intent(in ) :: mla + + ! local + integer n,nlevs,i,dm + + nlevs = mla%nlevel + dm = mla%dim + + do n=1,nlevs + do i=1,dm + call multifab_destroy(additive_sqrtLonsager_fc(n,i)) + end do + end do + + deallocate(additive_sqrtLonsager_fc) + + end subroutine destroy_additive_sqrtLonsager end module stochastic_mass_fluxdiv_module diff --git a/staggered_grid/src_reactDiff/main_driver.f90 b/staggered_grid/src_reactDiff/main_driver.f90 index a198a0fa..d6ecfd6a 100644 --- a/staggered_grid/src_reactDiff/main_driver.f90 +++ b/staggered_grid/src_reactDiff/main_driver.f90 @@ -358,6 +358,7 @@ subroutine main_driver() analyze_velocity=.false., & analyze_density=.true., & analyze_temperature=.false., & + analyze_Epot=.false., & structFactMultiplier = volume_factor) close(unit=un)