diff options
64 files changed, 1646 insertions, 177 deletions
diff --git a/.azure-pipelines.yml b/.azure-pipelines.yml index eee730896..dc71eea43 100644 --- a/.azure-pipelines.yml +++ b/.azure-pipelines.yml @@ -20,6 +20,10 @@ jobs: strategy: matrix: + cartesian1d: + WARPX_CI_REGULAR_CARTESIAN_1D: 'TRUE' + WARPX_CI_PSATD: 'FALSE' + WARPX_CI_OPENPMD: 'FALSE' cartesian2d: WARPX_CI_REGULAR_CARTESIAN_2D: 'TRUE' cartesian3d: diff --git a/.github/workflows/source/test_ci_matrix.sh b/.github/workflows/source/test_ci_matrix.sh index 104272eb4..cf9a69206 100755 --- a/.github/workflows/source/test_ci_matrix.sh +++ b/.github/workflows/source/test_ci_matrix.sh @@ -13,8 +13,10 @@ grep "\[" ci-tests.ini > ci_all_tests.txt export WARPX_CI_PSATD=TRUE # Concatenate the names of all elements in CI matrix into another test file -WARPX_CI_REGULAR_CARTESIAN_2D=TRUE python prepare_file_ci.py +WARPX_CI_REGULAR_CARTESIAN_1D=TRUE python prepare_file_ci.py grep "\[" ci-tests.ini > ci_matrix_elements.txt +WARPX_CI_REGULAR_CARTESIAN_2D=TRUE python prepare_file_ci.py +grep "\[" ci-tests.ini >> ci_matrix_elements.txt WARPX_CI_REGULAR_CARTESIAN_3D=TRUE python prepare_file_ci.py grep "\[" ci-tests.ini >> ci_matrix_elements.txt WARPX_CI_PYTHON_MAIN=TRUE python prepare_file_ci.py diff --git a/CMakeLists.txt b/CMakeLists.txt index d0dc02624..5c3574ed1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -60,8 +60,8 @@ option(WarpX_SENSEI "SENSEI in situ diagnostics" OFF) option(WarpX_QED "QED support (requires PICSAR)" ON) option(WarpX_QED_TABLE_GEN "QED table generation (requires PICSAR and Boost)" OFF) -set(WarpX_DIMS_VALUES 2 3 RZ) -set(WarpX_DIMS 3 CACHE STRING "Simulation dimensionality (2/3/RZ)") +set(WarpX_DIMS_VALUES 1 2 3 RZ) +set(WarpX_DIMS 3 CACHE STRING "Simulation dimensionality (1/2/3/RZ)") set_property(CACHE WarpX_DIMS PROPERTY STRINGS ${WarpX_DIMS_VALUES}) if(NOT WarpX_DIMS IN_LIST WarpX_DIMS_VALUES) message(FATAL_ERROR "WarpX_DIMS (${WarpX_DIMS}) must be one of ${WarpX_DIMS_VALUES}") @@ -263,6 +263,8 @@ if(WarpX_DIMS STREQUAL 3) target_compile_definitions(WarpX PUBLIC WARPX_DIM_3D) elseif(WarpX_DIMS STREQUAL 2) target_compile_definitions(WarpX PUBLIC WARPX_DIM_XZ) +elseif(WarpX_DIMS STREQUAL 1) + target_compile_definitions(WarpX PUBLIC WARPX_DIM_1D_Z) elseif(WarpX_DIMS STREQUAL RZ) target_compile_definitions(WarpX PUBLIC WARPX_DIM_RZ) endif() @@ -331,12 +333,10 @@ install(TARGETS ${WarpX_INSTALL_TARGET_NAMES} # simplified library alias # this is currently expected by Python bindings if(WarpX_LIB) - if(WarpX_DIMS STREQUAL 3) - set(lib_suffix "3d") - elseif(WarpX_DIMS STREQUAL 2) - set(lib_suffix "2d") - elseif(WarpX_DIMS STREQUAL RZ) + if(WarpX_DIMS STREQUAL RZ) set(lib_suffix "rz") + else() + set(lib_suffix "${WarpX_DIMS}d") endif() if(WIN32) set(mod_ext "dll") diff --git a/Docs/source/developers/dimensionality.rst b/Docs/source/developers/dimensionality.rst index ee98368d5..1ccdbba5d 100644 --- a/Docs/source/developers/dimensionality.rst +++ b/Docs/source/developers/dimensionality.rst @@ -13,6 +13,7 @@ Dimensions CMake Option ========== =================== **3D3V** ``WarpX_DIMS=3`` **2D3V** ``WarpX_DIMS=2`` +**1D3V** ``WarpX_DIMS=1`` **RZ** ``WarpX_DIMS=RZ`` ========== =================== @@ -23,31 +24,32 @@ Defines Depending on the build variant of WarpX, the following preprocessor macros will be set: -================== =========== =========== =========== -Macro 3D3V 2D3V RZ -================== =========== =========== =========== -``AMREX_SPACEDIM`` ``3`` ``2`` ``2`` -``WARPX_DIM_3D`` **defined** *undefined* *undefined* -``WARPX_DIM_XZ`` *undefined* **defined** *undefined* -``WARPX_DIM_RZ`` *undefined* *undefined* **defined** -================== =========== =========== =========== +================== =========== =========== =========== =========== +Macro 3D3V 2D3V 1D3V RZ +================== =========== =========== =========== =========== +``AMREX_SPACEDIM`` ``3`` ``2`` ``1`` ``2`` +``WARPX_DIM_3D`` **defined** *undefined* *undefined* *undefined* +``WARPX_DIM_1D_Z`` *undefined* *undefined* **defined** *undefined* +``WARPX_DIM_XZ`` *undefined* **defined** *undefined* *undefined* +``WARPX_DIM_RZ`` *undefined* *undefined* *undefined* **defined** +================== =========== =========== =========== =========== At the same time, the following conventions will apply: -==================== =========== =========== =========== -**Convention** **3D3V** **2D3V** **RZ** --------------------- ----------- ----------- ----------- +==================== =========== =========== =========== =========== +**Convention** **3D3V** **2D3V** **1D3V** **RZ** +-------------------- ----------- ----------- ----------- ----------- *Fields* ------------------------------------------------------------ -AMReX Box dimensions ``3`` ``2`` ``2`` -WarpX axis labels ``x, y, z`` ``x, z`` ``x, z`` --------------------- ----------- ----------- ----------- +------------------------------------------------------------------------ +AMReX Box dimensions ``3`` ``2`` ``1`` ``2`` +WarpX axis labels ``x, y, z`` ``x, z`` ``z`` ``x, z`` +-------------------- ----------- ----------- ----------- ----------- *Particles* ------------------------------------------------------------ -AMReX AoS ``.pos()`` ``0, 1, 2`` ``0, 1`` ``0, 1`` -WarpX position names ``x, y, z`` ``x, z`` ``r, z`` -extra SoA attribute ``theta`` -==================== =========== =========== =========== +------------------------------------------------------------------------ +AMReX AoS ``.pos()`` ``0, 1, 2`` ``0, 1`` ``0`` ``0, 1`` +WarpX position names ``x, y, z`` ``x, z`` ``z`` ``r, z`` +extra SoA attribute ``theta`` +==================== =========== =========== =========== =========== Please see the following sections for particle AoS and SoA details. @@ -55,3 +57,4 @@ Conventions ----------- In 2D3V, we assume that the position of a particle in ``y`` is equal to ``0``. +In 1D3V, we assume that the position of a particle in ``x`` and ``y`` is equal to ``0``. diff --git a/Examples/Modules/laser_injection/analysis_1d.py b/Examples/Modules/laser_injection/analysis_1d.py new file mode 100755 index 000000000..6606938df --- /dev/null +++ b/Examples/Modules/laser_injection/analysis_1d.py @@ -0,0 +1,175 @@ +#! /usr/bin/env python + +# Copyright 2021 Prabhat Kumar, Remi Lehe +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +# This file is part of the WarpX automated test suite. Its purpose is to test the +# injection of a Gaussian laser pulse from an antenna in a 1D simulation. +# The test calculates the envelope of each component of the laser pulse at the end of +# the simulation and it compares it with theory. It also checks that the +# central frequency of the Fourier transform is the expected one. + +import yt +import sys +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import numpy as np +from scipy.signal import hilbert +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# Maximum acceptable error for this test +relative_error_threshold = 0.05 + +# A small number +small_num = 1.0e-8 + +# Physical parameters +um = 1.e-6 +fs = 1.e-15 +c = 299792458 + +# Parameters of the gaussian beam +wavelength = 1.*um +w0 = 5.*um +tt = 10.*fs +t_c = 24.*fs +E_max = 4e12 + +# laser direction +dir_vector = np.array([0,0,1.0]) +dir_vector /= np.linalg.norm(dir_vector) + + +# polarization vector +pol_vector = np.array([1.0,1.0,0.0]) +pol_vector /= np.linalg.norm(pol_vector) + +# Calculates the envelope of a Gaussian beam +def gauss_env(T,Z): + '''Function to compute the theory for the envelope + ''' + inv_tau2 = 1./tt/tt + exp_arg = - inv_tau2 / c/c * (Z-T*c)*(Z-T*c) + return E_max * np.real(np.exp(exp_arg)) + +# Checks envelope and central frequency for a given laser component +def check_component(data, component, t_env_theory, coeff,Z,dz): + print("*** Checking " + component + " ***") + field = data['boxlib', component].v.squeeze() + env = abs(hilbert(field)) + + env_theory = t_env_theory*np.abs(coeff) + + # Plot results + fig = plt.figure(figsize=(12,6)) + + ax1 = fig.add_subplot(221) + ax1.set_title('PIC field') + ax1.plot(Z,field) + + ax2 = fig.add_subplot(222) + ax2.set_title('PIC envelope') + ax2.plot(Z,env) + + ax3 = fig.add_subplot(223) + ax3.set_title('Theory envelope') + ax3.plot(Z,env_theory, label="theory") + ax3.plot(Z,env, label="simulation") + ax3.legend(loc="upper right") + + ax4 = fig.add_subplot(224) + ax4.set_title('Difference') + ax4.plot(Z,env-env_theory) + + plt.tight_layout() + plt.savefig("plt_" + component + ".png", bbox_inches='tight') + + if(np.abs(coeff) < small_num): + is_field_zero = np.sum(np.abs(env)) < small_num + if is_field_zero : + print("[OK] Field component expected to be 0 is ~ 0") + else : + print("[FAIL] Field component expected to be 0 is NOT ~ 0") + assert(is_field_zero) + print("******\n") + return + + fft_field = np.fft.fft(field) + + freq_cols = np.fft.fftfreq(fft_field.shape[0],dz/c) + + pos_max = np.unravel_index(np.abs(fft_field).argmax(), fft_field.shape) + + freq = np.abs(freq_cols[pos_max[0]]) + exp_freq = c/wavelength + + relative_error_freq = np.abs(freq-exp_freq)/exp_freq + is_freq_ok = relative_error_freq < relative_error_threshold + if is_freq_ok : + print("[OK] Relative error frequency: {:6.3f} %".format(relative_error_freq*100)) + else : + print("[FAIL] Relative error frequency: {:6.3f} %".format(relative_error_freq*100)) + assert(is_freq_ok) + + print("******\n") + + relative_error_env = np.sum(np.abs(env-env_theory)) / np.sum(np.abs(env_theory)) + is_env_ok = relative_error_env < relative_error_threshold + if is_env_ok : + print("[OK] Relative error envelope: {:6.3f} %".format(relative_error_env*100)) + else : + print("[FAIL] Relative error envelope: {:6.3f} %".format(relative_error_env*100)) + assert(is_env_ok) + +def check_laser(filename): + ds = yt.load(filename) + + # yt 4.0+ has rounding issues with our domain data: + # RuntimeError: yt attempted to read outside the boundaries + # of a non-periodic domain along dimension 0. + if 'force_periodicity' in dir(ds): ds.force_periodicity() + + z = np.linspace( + ds.domain_left_edge[0].v, + ds.domain_right_edge[0].v, + ds.domain_dimensions[0]) + + dz = (ds.domain_right_edge[0].v-ds.domain_left_edge[0].v)/(ds.domain_dimensions[0]-1) + + # Compute the theory for envelope + env_theory = gauss_env(+t_c-ds.current_time.to_value(),z)+gauss_env(-t_c+ds.current_time.to_value(),z) + + # Read laser field in PIC simulation, and compute envelope + all_data_level_0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) + + b_vector = np.cross(dir_vector, pol_vector) + + components = ["Ex", "Ey", "Ez", "Bx", "By", "Bz"] + coeffs = [ + pol_vector[0], + pol_vector[1], + pol_vector[2], + b_vector[0], + b_vector[1], + b_vector[2]] + + field_facts = [1, 1, 1, 1/c, 1/c, 1/c] + + for comp, coeff, field_fact in zip(components, coeffs, field_facts): + check_component(all_data_level_0, comp, field_fact*env_theory, coeff, z, dz) + +def main(): + filename_end = sys.argv[1] + + check_laser(filename_end) + + test_name = filename_end[:-9] # Could also be os.path.split(os.getcwd())[1] + checksumAPI.evaluate_checksum(test_name, filename_end) + +if __name__ == "__main__": + main() diff --git a/Examples/Modules/laser_injection/inputs_1d_rt b/Examples/Modules/laser_injection/inputs_1d_rt new file mode 100644 index 000000000..eb947039c --- /dev/null +++ b/Examples/Modules/laser_injection/inputs_1d_rt @@ -0,0 +1,62 @@ +# Maximum number of time steps +max_step = 240 + +# number of grid points +amr.n_cell = 352 + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +amr.max_grid_size = 32 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -15.e-6 # physical domain +geometry.prob_hi = 15.e-6 + +boundary.field_lo = pec +boundary.field_hi = pec + +warpx.serialize_ics = 1 + +# Verbosity +warpx.verbose = 1 + +# Algorithms +algo.current_deposition = esirkepov +warpx.use_filter = 0 + +# CFL +warpx.cfl = 0.9 + +# Order of particle shape factors +algo.particle_shape = 1 + +# Laser +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0.e-6 0.e-6 0.e-6 # This point is on the laser plane +laser1.direction = 0. 0. 1. # The plane normal direction +laser1.polarization = 1. 1. 0. # The main polarization vector +laser1.e_max = 4.e12 # Maximum amplitude of the laser field (in V/m) +laser1.wavelength = 1.0e-6 # The wavelength of the laser (in meters) +laser1.profile_waist = 5.e-6 # The waist of the laser (in meters) +laser1.profile_duration = 10.e-15 # The duration of the laser (in seconds) +laser1.profile_t_peak = 24.e-15 # The time at which the laser reaches its peak (in seconds) +laser1.profile_focal_distance = 13.109e-6 # Focal distance from the antenna (in meters) + # With this focal distance the laser is at focus + # at the end of the simulation. + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 20 +diag1.diag_type = Full + +# Moving window +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 # in units of the speed of light +warpx.start_moving_window_step = 20 +warpx.end_moving_window_step = 200 diff --git a/Examples/Tests/Langmuir/analysis_langmuir_multi_1d.py b/Examples/Tests/Langmuir/analysis_langmuir_multi_1d.py new file mode 100755 index 000000000..5cb52de23 --- /dev/null +++ b/Examples/Tests/Langmuir/analysis_langmuir_multi_1d.py @@ -0,0 +1,114 @@ +#! /usr/bin/env python + +# Copyright 2019-2021 Jean-Luc Vay, Maxence Thevenet, Remi Lehe, Prabhat Kumar +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + + +# This is a script that analyses the simulation results from +# the script `inputs.multi.rt`. This simulates a 1D periodic plasma wave. +# The electric field in the simulation is given (in theory) by: +# $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\sin(k_z z)\sin( \omega_p t)$$ +import sys +import re +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import yt +yt.funcs.mylog.setLevel(50) +import numpy as np +from scipy.constants import e, m_e, epsilon_0, c +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# this will be the name of the plot file +fn = sys.argv[1] + +# Parse test name and check if current correction (psatd.current_correction=1) is applied +current_correction = True if re.search( 'current_correction', fn ) else False + +# Parse test name and check if Vay current deposition (algo.current_deposition=vay) is used +vay_deposition = True if re.search( 'Vay_deposition', fn ) else False + +# Parameters (these parameters must match the parameters in `inputs.multi.rt`) +epsilon = 0.01 +n = 4.e24 +n_osc_z = 2 +zmin = -20e-6; zmax = 20.e-6; Nz = 128 + +# Wave vector of the wave +kz = 2.*np.pi*n_osc_z/(zmax-zmin) +# Plasma frequency +wp = np.sqrt((n*e**2)/(m_e*epsilon_0)) + +k = {'Ez':kz} +cos = {'Ez':(1,1,0)} + +def get_contribution( is_cos, k ): + du = (zmax-zmin)/Nz + u = zmin + du*( 0.5 + np.arange(Nz) ) + if is_cos == 1: + return( np.cos(k*u) ) + else: + return( np.sin(k*u) ) + +def get_theoretical_field( field, t ): + amplitude = epsilon * (m_e*c**2*k[field])/e * np.sin(wp*t) + cos_flag = cos[field] + z_contribution = get_contribution( cos_flag[2], kz ) + + E = amplitude * z_contribution + + return( E ) + +# Read the file +ds = yt.load(fn) +t0 = ds.current_time.to_value() +data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) +# Check the validity of the fields +error_rel = 0 +for field in ['Ez']: + E_sim = data[field].to_ndarray()[:,0,0] + E_th = get_theoretical_field(field, t0) + max_error = abs(E_sim-E_th).max()/abs(E_th).max() + print('%s: Max error: %.2e' %(field,max_error)) + error_rel = max( error_rel, max_error ) + +# Plot the last field from the loop (Ez at iteration 80) +plt.subplot2grid( (1,2), (0,0) ) +plt.plot( E_sim ) +#plt.colorbar() +plt.title('Ez, last iteration\n(simulation)') +plt.subplot2grid( (1,2), (0,1) ) +plt.plot( E_th ) +#plt.colorbar() +plt.title('Ez, last iteration\n(theory)') +plt.tight_layout() +plt.savefig('langmuir_multi_1d_analysis.png') + +tolerance_rel = 0.05 + +print("error_rel : " + str(error_rel)) +print("tolerance_rel: " + str(tolerance_rel)) + +assert( error_rel < tolerance_rel ) + +# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when +# current correction (psatd.do_current_correction=1) is applied or when +# Vay current deposition (algo.current_deposition=vay) is used +if current_correction or vay_deposition: + rho = data['rho' ].to_ndarray() + divE = data['divE'].to_ndarray() + error_rel = np.amax( np.abs( divE - rho/epsilon_0 ) ) / np.amax( np.abs( rho/epsilon_0 ) ) + tolerance = 1.e-9 + print("Check charge conservation:") + print("error_rel = {}".format(error_rel)) + print("tolerance = {}".format(tolerance)) + assert( error_rel < tolerance ) + +test_name = fn[:-9] # Could also be os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/Langmuir/inputs_1d_multi_rt b/Examples/Tests/Langmuir/inputs_1d_multi_rt new file mode 100644 index 000000000..9ed38d790 --- /dev/null +++ b/Examples/Tests/Langmuir/inputs_1d_multi_rt @@ -0,0 +1,81 @@ +# Maximum number of time steps +max_step = 80 + +# number of grid points +amr.n_cell = 128 + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +amr.max_grid_size = 64 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.prob_lo = -20.e-6 # physical domain +geometry.prob_hi = 20.e-6 + +# Boundary condition +boundary.field_lo = periodic +boundary.field_hi = periodic + +warpx.serialize_ics = 1 + +# Verbosity +warpx.verbose = 1 + +# Algorithms +algo.field_gathering = energy-conserving +warpx.use_filter = 0 + +# Order of particle shape factors +algo.particle_shape = 1 + +# CFL +warpx.cfl = 0.8 + +# Parameters for the plasma wave +my_constants.epsilon = 0.01 +my_constants.n0 = 2.e24 # electron and positron densities, #/m^3 +my_constants.wp = sqrt(2.*n0*q_e**2/(epsilon0*m_e)) # plasma frequency +my_constants.kp = wp/clight # plasma wavenumber +my_constants.k = 2.*pi/20.e-6 # perturbation wavenumber +# Note: kp is calculated in SI for a density of 4e24 (i.e. 2e24 electrons + 2e24 positrons) +# k is calculated so as to have 2 periods within the 40e-6 wide box. + +# Particles +particles.species_names = electrons positrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = 2 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.profile = constant +electrons.density = n0 # number of electrons per m^3 +electrons.momentum_distribution_type = parse_momentum_function +electrons.momentum_function_ux(x,y,z) = "epsilon * k/kp * sin(k*x) * cos(k*y) * cos(k*z)" +electrons.momentum_function_uy(x,y,z) = "epsilon * k/kp * cos(k*x) * sin(k*y) * cos(k*z)" +electrons.momentum_function_uz(x,y,z) = "epsilon * k/kp * cos(k*x) * cos(k*y) * sin(k*z)" + +positrons.charge = q_e +positrons.mass = m_e +positrons.injection_style = "NUniformPerCell" +positrons.num_particles_per_cell_each_dim = 2 +positrons.zmin = -20.e-6 +positrons.zmax = 20.e-6 + +positrons.profile = constant +positrons.density = n0 # number of positrons per m^3 +positrons.momentum_distribution_type = parse_momentum_function +positrons.momentum_function_ux(x,y,z) = "-epsilon * k/kp * sin(k*x) * cos(k*y) * cos(k*z)" +positrons.momentum_function_uy(x,y,z) = "-epsilon * k/kp * cos(k*x) * sin(k*y) * cos(k*z)" +positrons.momentum_function_uz(x,y,z) = "-epsilon * k/kp * cos(k*x) * cos(k*y) * sin(k*z)" + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 40 +diag1.diag_type = Full diff --git a/GNUmakefile b/GNUmakefile index 0842c0696..7ddded504 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -8,7 +8,8 @@ DEBUG = FALSE WARN_ALL = TRUE #WARN_ERROR=TRUE -#DIM = 2 +#DIM = 1 +#DIM = 2 DIM = 3 QED = TRUE diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_1d.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_1d.json new file mode 100644 index 000000000..974b7cdb5 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_1d.json @@ -0,0 +1,31 @@ +{ + "electrons": { + "particle_cpu": 0.0, + "particle_id": 0.0, + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 0.0, + "particle_position_x": 0.0, + "particle_weight": 0.0 + }, + "lev=0": { + "Bx": 0.0, + "By": 0.0, + "Bz": 0.0, + "Ex": 0.0, + "Ey": 0.0, + "Ez": 124191078376.51929, + "jx": 0.0, + "jy": 0.0, + "jz": 47472967251388.14 + }, + "positrons": { + "particle_cpu": 0.0, + "particle_id": 0.0, + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 0.0, + "particle_position_x": 0.0, + "particle_weight": 0.0 + } +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/LaserInjection_1d.json b/Regression/Checksum/benchmarks_json/LaserInjection_1d.json new file mode 100644 index 000000000..e40fd95ba --- /dev/null +++ b/Regression/Checksum/benchmarks_json/LaserInjection_1d.json @@ -0,0 +1,13 @@ +{ + "lev=0": { + "Bx": 374663.2330698118, + "By": 374800.98321004683, + "Bz": 0.0, + "Ex": 111565176079461.14, + "Ey": 111523937757159.5, + "Ez": 0.0, + "jx": 116159381829.31923, + "jy": 116159381829.31923, + "jz": 0.0 + } +}
\ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 30cc290dd..07540bb25 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -687,6 +687,25 @@ analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py analysisOutputImage = langmuir_multi_2d_analysis.png tolerance = 1.e-14 +[Langmuir_multi_1d] +buildDir = . +inputFile = Examples/Tests/Langmuir/inputs_1d_multi_rt +runtime_params = algo.current_deposition=esirkepov diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz +dim = 1 +addToCompileString = QED=FALSE +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons positrons +analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_1d.py +analysisOutputImage = langmuir_multi_1d_analysis.png +tolerance = 1.e-14 + [Langmuir_multi_rz] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rz_rt @@ -820,6 +839,23 @@ compareParticles = 0 analysisRoutine = Examples/Modules/laser_injection/analysis_2d.py tolerance = 1.e-14 +[LaserInjection_1d] +buildDir = . +inputFile = Examples/Modules/laser_injection/inputs_1d_rt +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_ics=1 +dim = 1 +addToCompileString = QED=FALSE +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Modules/laser_injection/analysis_1d.py +tolerance = 1.e-14 + [LaserAcceleration] buildDir = . inputFile = Examples/Physics_applications/laser_acceleration/inputs_3d diff --git a/Regression/prepare_file_ci.py b/Regression/prepare_file_ci.py index bb94cb992..46c9bf8b1 100644 --- a/Regression/prepare_file_ci.py +++ b/Regression/prepare_file_ci.py @@ -15,6 +15,7 @@ import os # Get relevant environment variables arch = os.environ.get('WARPX_TEST_ARCH', 'CPU') +ci_regular_cartesian_1d = os.environ.get('WARPX_CI_REGULAR_CARTESIAN_1D') == 'TRUE' ci_regular_cartesian_2d = os.environ.get('WARPX_CI_REGULAR_CARTESIAN_2D') == 'TRUE' ci_regular_cartesian_3d = os.environ.get('WARPX_CI_REGULAR_CARTESIAN_3D') == 'TRUE' ci_psatd = os.environ.get('WARPX_CI_PSATD', 'TRUE') == 'TRUE' @@ -109,6 +110,15 @@ def select_tests(blocks, match_string_list, do_test): blocks = [ block for block in blocks if match_string in block ] return blocks +if ci_regular_cartesian_1d: + test_blocks = select_tests(test_blocks, ['dim = 1'], True) + test_blocks = select_tests(test_blocks, ['USE_RZ=TRUE'], False) + test_blocks = select_tests(test_blocks, ['PYTHON_MAIN=TRUE'], False) + test_blocks = select_tests(test_blocks, ['PRECISION=FLOAT', 'USE_SINGLE_PRECISION_PARTICLES=TRUE'], False) + test_blocks = select_tests(test_blocks, ['useMPI = 0'], False) + test_blocks = select_tests(test_blocks, ['QED=TRUE'], False) + test_blocks = select_tests(test_blocks, ['USE_EB=TRUE'], False) + if ci_regular_cartesian_2d: test_blocks = select_tests(test_blocks, ['dim = 2'], True) test_blocks = select_tests(test_blocks, ['USE_RZ=TRUE'], False) @@ -119,8 +129,7 @@ if ci_regular_cartesian_2d: test_blocks = select_tests(test_blocks, ['USE_EB=TRUE'], False) if ci_regular_cartesian_3d: - test_blocks = select_tests(test_blocks, ['dim = 2'], False) - test_blocks = select_tests(test_blocks, ['USE_RZ=TRUE'], False) + test_blocks = select_tests(test_blocks, ['dim = 3'], True) test_blocks = select_tests(test_blocks, ['PYTHON_MAIN=TRUE'], False) test_blocks = select_tests(test_blocks, ['PRECISION=FLOAT', 'USE_SINGLE_PRECISION_PARTICLES=TRUE'], False) test_blocks = select_tests(test_blocks, ['useMPI = 0'], False) diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 441acb743..47e4cce62 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -184,7 +184,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { +#if (AMREX_SPACEDIM >= 2) int jdim = (idim+1) % AMREX_SPACEDIM; +#endif #if (AMREX_SPACEDIM == 3) int kdim = (idim+2) % AMREX_SPACEDIM; #endif @@ -198,6 +200,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n { direct_faces.push_back(kv.first); } +#if (AMREX_SPACEDIM >= 2) else if (amrex::grow(grid_box, jdim, ncell).intersects(box)) { side_faces.push_back(kv.first); @@ -227,8 +230,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n { corners.push_back(kv.first); } +#endif } +#if (AMREX_SPACEDIM >= 2) for (auto gid : corners) { const Box& grid_box = grids[gid]; @@ -261,6 +266,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n amrex::Abort("SigmaBox::SigmaBox(): corners, how did this happen?\n"); } } +#endif #if (AMREX_SPACEDIM == 3) for (auto gid : side_side_edges) @@ -302,6 +308,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n } #endif +#if (AMREX_SPACEDIM >= 2) for (auto gid : side_faces) { const Box& grid_box = grids[gid]; @@ -317,6 +324,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); } } +#endif for (auto gid : direct_faces) { @@ -369,7 +377,12 @@ SigmaBox::ComputePMLFactorsB (const Real* a_dx, Real dt) N[idim] = sigma_star[idim].size(); dx[idim] = a_dx[idim]; } - amrex::ParallelFor(amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), + amrex::ParallelFor( +#if (AMREX_SPACEDIM >= 2) + amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), +#else + N[0], +#endif [=] AMREX_GPU_DEVICE (int i) noexcept { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -398,7 +411,12 @@ SigmaBox::ComputePMLFactorsE (const Real* a_dx, Real dt) N[idim] = sigma[idim].size(); dx[idim] = a_dx[idim]; } - amrex::ParallelFor(amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), + amrex::ParallelFor( +#if (AMREX_SPACEDIM >= 2) + amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), +#else + N[0], +#endif [=] AMREX_GPU_DEVICE (int i) noexcept { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -524,6 +542,8 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g IntVect ngFFT = IntVect(ngFFt_x, ngFFt_y, ngFFt_z); #elif (AMREX_SPACEDIM == 2) IntVect ngFFT = IntVect(ngFFt_x, ngFFt_z); +#elif (AMREX_SPACEDIM == 1) + IntVect ngFFT = IntVect(ngFFt_z); #endif // Set the number of guard cells to the maximum of each field diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 9c7435b6c..4458f1e2a 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -80,6 +80,9 @@ PEC::ApplyPECtoEfield (std::array<amrex::MultiFab*, 3> Efield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 0; PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -90,6 +93,9 @@ PEC::ApplyPECtoEfield (std::array<amrex::MultiFab*, 3> Efield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 1; PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -100,6 +106,9 @@ PEC::ApplyPECtoEfield (std::array<amrex::MultiFab*, 3> Efield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 2; PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -163,6 +172,9 @@ PEC::ApplyPECtoBfield (std::array<amrex::MultiFab*, 3> Bfield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 0; PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -173,6 +185,9 @@ PEC::ApplyPECtoBfield (std::array<amrex::MultiFab*, 3> Bfield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 1; PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -183,6 +198,9 @@ PEC::ApplyPECtoBfield (std::array<amrex::MultiFab*, 3> Bfield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 2; PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, diff --git a/Source/BoundaryConditions/WarpX_PML_kernels.H b/Source/BoundaryConditions/WarpX_PML_kernels.H index ca6b594ac..c434e09f3 100644 --- a/Source/BoundaryConditions/WarpX_PML_kernels.H +++ b/Source/BoundaryConditions/WarpX_PML_kernels.H @@ -27,6 +27,12 @@ void warpx_damp_pml_ex (int i, int j, int k, Array4<Real> const& Ex, int xlo, int ylo, int zlo, const bool dive_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Ex, Ex_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -98,6 +104,12 @@ void warpx_damp_pml_ey (int i, int j, int k, Array4<Real> const& Ey, int xlo, int ylo, int zlo, const bool dive_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Ey, Ey_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, dive_cleaning); @@ -166,6 +178,12 @@ void warpx_damp_pml_ez (int i, int j, int k, Array4<Real> const& Ez, int xlo, int ylo, int zlo, const bool dive_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Ez, Ez_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -237,6 +255,12 @@ void warpx_damp_pml_bx (int i, int j, int k, Array4<Real> const& Bx, int xlo, int ylo, int zlo, const bool divb_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Bx, Bx_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -308,6 +332,12 @@ void warpx_damp_pml_by (int i, int j, int k, Array4<Real> const& By, int xlo, int ylo, int zlo, const bool divb_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, By, By_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, divb_cleaning); @@ -376,6 +406,12 @@ void warpx_damp_pml_bz (int i, int j, int k, Array4<Real> const& Bz, int xlo, int ylo, int zlo, const bool divb_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Bz, Bz_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -446,6 +482,12 @@ void warpx_damp_pml_scalar (int i, int j, int k, Array4<Real> const& arr, const Real* const sigma_star_fac_z, int xlo, int ylo, int zlo) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, arr, arr_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index 2f2231c0c..236516f07 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -301,6 +301,7 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) / dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]) ) ); // Take the max of 0 and num_zcells_lab int Nz_lab = std::max( 0, num_zcells_lab ); +#if (AMREX_SPACEDIM >= 2) // Number of lab-frame cells in x-direction at level, lev const int num_xcells_lab = static_cast<int>( floor ( ( diag_dom.hi(0) - diag_dom.lo(0) ) @@ -308,6 +309,7 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) ) ); // Take the max of 0 and num_ycells_lab int Nx_lab = std::max( 0, num_xcells_lab); +#endif #if (AMREX_SPACEDIM == 3) // Number of lab-frame cells in the y-direction at level, lev const int num_ycells_lab = static_cast<int>( floor ( @@ -317,8 +319,10 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) // Take the max of 0 and num_xcells_lab int Ny_lab = std::max( 0, num_ycells_lab ); m_snapshot_ncells_lab[i_buffer] = {Nx_lab, Ny_lab, Nz_lab}; -#else +#elif (AMREX_SPACEDIM == 2) m_snapshot_ncells_lab[i_buffer] = {Nx_lab, Nz_lab}; +#else + m_snapshot_ncells_lab[i_buffer] = amrex::IntVect({Nz_lab}); #endif } diff --git a/Source/Diagnostics/BackTransformedDiagnostic.cpp b/Source/Diagnostics/BackTransformedDiagnostic.cpp index 17b59d3ae..c2d255d06 100644 --- a/Source/Diagnostics/BackTransformedDiagnostic.cpp +++ b/Source/Diagnostics/BackTransformedDiagnostic.cpp @@ -160,8 +160,10 @@ namespace // Create a 3D, nx x ny x nz dataspace. #if (AMREX_SPACEDIM == 3) hsize_t dims[3] = {nx, ny, nz}; -#else +#elif (AMREX_SPACEDIM == 2) hsize_t dims[3] = {nx, nz}; +#else + hsize_t dims[3] = {nz}; #endif hid_t grid_space = H5Screate_simple(AMREX_SPACEDIM, dims, NULL); @@ -422,11 +424,15 @@ namespace shift[0] = lo_x; shift[1] = lo_y; shift[2] = lo_z; -#else +#elif (AMREX_SPACEDIM == 2) hsize_t slab_offsets[2], slab_dims[2]; int shift[2]; shift[0] = lo_x; shift[1] = lo_z; +#else + hsize_t slab_offsets[1], slab_dims[1]; + int shift[1]; + shift[0] = lo_z; #endif hid_t slab_dataspace; @@ -584,13 +590,19 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, m_dz_lab_ = PhysConst::c * m_dt_boost_ * m_inv_beta_boost_ * m_inv_gamma_boost_; m_inv_dz_lab_ = 1.0_rt / m_dz_lab_; int Nz_lab = static_cast<unsigned>((zmax_lab - zmin_lab) * m_inv_dz_lab_); +#if (AMREX_SPACEDIM >= 2) int Nx_lab = geom.Domain().length(0); +#endif #if (AMREX_SPACEDIM == 3) int Ny_lab = geom.Domain().length(1); IntVect prob_ncells_lab = {Nx_lab, Ny_lab, Nz_lab}; -#else +#elif (AMREX_SPACEDIM == 2) // Ny_lab = 1; IntVect prob_ncells_lab = {Nx_lab, Nz_lab}; +#else + // Nx_lab = 1; + // Ny_lab = 1; + IntVect prob_ncells_lab({Nz_lab}); #endif writeMetaData(); @@ -640,7 +652,6 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, for (int i = 0; i < N_slice_snapshots; ++i) { - IntVect slice_ncells_lab ; // To construct LabFrameSlice(), the location of lo() and hi() of the // reduced diag is computed using the user-defined values of the @@ -657,6 +668,7 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, ( (1._rt+m_beta_boost_)*m_gamma_boost_); auto Nz_slice_lab = static_cast<int>( (zmax_slice_lab - zmin_slice_lab) * m_inv_dz_lab_); +#if (AMREX_SPACEDIM >= 2) auto Nx_slice_lab = static_cast<int>( (current_slice_hi[0] - current_slice_lo[0] ) / geom.CellSize(0)); @@ -664,6 +676,7 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, // if the x-dimension is reduced, increase total_cells by 1 // to be consistent with the number of cells created for the output. if (Nx_lab != Nx_slice_lab) Nx_slice_lab++; +#endif #if (AMREX_SPACEDIM == 3) auto Ny_slice_lab = static_cast<int>( (current_slice_hi[1] - current_slice_lo[1]) / @@ -672,9 +685,11 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, // if the y-dimension is reduced, increase total_cells by 1 // to be consistent with the number of cells created for the output. if (Ny_lab != Ny_slice_lab) Ny_slice_lab++; - slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab}; + amrex::IntVect slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab}; +#elif (AMREX_SPACEDIM == 2) + amrex::IntVect slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab}; #else - slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab}; + amrex::IntVect slice_ncells_lab({Nz_slice_lab}); #endif IntVect slice_lo(AMREX_D_DECL(0,0,0)); @@ -1213,7 +1228,11 @@ createLabFrameDirectories() { const auto lo = lbound(m_buff_box_); for (int comp = 0; comp < m_ncomp_to_dump_; ++comp) { output_create_field(m_file_name, m_mesh_field_names_[comp], +#if ( AMREX_SPACEDIM >= 2 ) m_prob_ncells_lab_[0], +#else + 1, +#endif #if ( AMREX_SPACEDIM == 3 ) m_prob_ncells_lab_[1], #else @@ -1387,8 +1406,10 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, const int icomp = field_map_ptr[n]; #if (AMREX_SPACEDIM == 3) buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp); -#else +#elif (AMREX_SPACEDIM == 2) buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp); +#else + buf_arr(k_lab,j,k,n) = tmp_arr(i,j,k,icomp); #endif } ); @@ -1424,8 +1445,10 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, const int icomp = field_map_ptr[n]; #if (AMREX_SPACEDIM == 3) buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp); -#else +#elif (AMREX_SPACEDIM == 2) buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp); +#else + buf_arr(k_lab,j,k,n) = tmp_arr(i,j,k,icomp); #endif }); } @@ -1531,8 +1554,10 @@ AddPartDataToParticleBuffer( int* const AMREX_RESTRICT IndexLocation = IndexForPartCopy.dataPtr(); // Compute extent of the reduced domain +/- user-defined physical width +#if (AMREX_SPACEDIM >= 2) Real const xmin = m_diag_domain_lab_.lo(0)-m_particle_slice_dx_lab_; Real const xmax = m_diag_domain_lab_.hi(0)+m_particle_slice_dx_lab_; +#endif #if (AMREX_SPACEDIM == 3) Real const ymin = m_diag_domain_lab_.lo(1)-m_particle_slice_dx_lab_; Real const ymax = m_diag_domain_lab_.hi(1)+m_particle_slice_dx_lab_; @@ -1544,8 +1569,11 @@ AddPartDataToParticleBuffer( [=] AMREX_GPU_DEVICE(int i) { Flag[i] = 0; +#if (AMREX_SPACEDIM >= 2) if ( x_temp[i] >= (xmin) && - x_temp[i] <= (xmax) ) { + x_temp[i] <= (xmax) ) +#endif + { #if (AMREX_SPACEDIM == 3) if (y_temp[i] >= (ymin) && y_temp[i] <= (ymax) ) diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index 9df599fa7..6fef39ce8 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -109,8 +109,10 @@ Diagnostics::BaseReadParameters () if (warpx.do_moving_window) { #if (AMREX_SPACEDIM == 3) amrex::Vector<int> dim_map {0, 1, 2}; -#else +#elif (AMREX_SPACEDIM == 2) amrex::Vector<int> dim_map {0, 2}; +#else + amrex::Vector<int> dim_map {2}; #endif if (warpx.boost_direction[ dim_map[warpx.moving_window_dir] ] == 1) { // Convert user-defined lo and hi for diagnostics to account for boosted-frame diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index 3490c4755..4bcf233fb 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -496,6 +496,11 @@ FullDiagnostics::MovingWindowAndGalileanDomainShift (int step) new_lo[1] = current_lo[1] + warpx.m_galilean_shift[2]; new_hi[1] = current_hi[1] + warpx.m_galilean_shift[2]; } +#elif (AMREX_SPACEDIM == 1 ) + { + new_lo[0] = current_lo[0] + warpx.m_galilean_shift[2]; + new_hi[0] = current_hi[0] + warpx.m_galilean_shift[2]; + } #endif // Update RealBox of geometry with galilean-shifted boundary for (int lev = 0; lev < nmax_lev; ++lev) { diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp index 4896176c2..0c2cf7e95 100644 --- a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp @@ -60,6 +60,16 @@ BeamRelevant::BeamRelevant (std::string rd_name) // 12,13: emittance x,z // 14: charge m_data.resize(15, 0.0_rt); +#elif (defined WARPX_DIM_1D_Z) + // 0 : mean z + // 1,2,3 : mean px,py,pz + // 4 : gamma + // 5 : rms z + // 6,7,8 : rms px,py,pz + // 9 : rms gamma + // 10 : emittance z + // 11 : charge + m_data.resize(12, 0.0_rt); #endif if (ParallelDescriptor::IOProcessor()) @@ -112,6 +122,23 @@ BeamRelevant::BeamRelevant (std::string rd_name) ofs << "[" << c++ << "]emittance_x(m)"; ofs << m_sep; ofs << "[" << c++ << "]emittance_z(m)"; ofs << m_sep; ofs << "[" << c++ << "]charge(C)"; ofs << std::endl; +#elif (defined WARPX_DIM_1D_Z) + int c = 0; + ofs << "#"; + ofs << "[" << c++ << "]step()"; ofs << m_sep; + ofs << "[" << c++ << "]time(s)"; ofs << m_sep; + ofs << "[" << c++ << "]z_mean(m)"; ofs << m_sep; + ofs << "[" << c++ << "]px_mean(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]py_mean(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]pz_mean(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]gamma_mean()"; ofs << m_sep; + ofs << "[" << c++ << "]z_rms(m)"; ofs << m_sep; + ofs << "[" << c++ << "]px_rms(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]py_rms(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]pz_rms(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]gamma_rms()"; ofs << m_sep; + ofs << "[" << c++ << "]emittance_z(m)"; ofs << m_sep; + ofs << "[" << c++ << "]charge(C)"; ofs << std::endl; #endif // close file ofs.close(); @@ -143,6 +170,8 @@ void BeamRelevant::ComputeDiags (int step) int const index_z = 2; #elif (defined WARPX_DIM_XZ || defined WARPX_DIM_RZ) int const index_z = 1; +#elif (defined WARPX_DIM_1D_Z) + int const index_z = 0; #endif // loop over species @@ -250,10 +279,16 @@ void BeamRelevant::ComputeDiags (int step) const ParticleReal p_pos0 = p.pos(0); const ParticleReal p_w = p.rdata(PIdx::w); -#if (defined WARPX_DIM_RZ) +#if (defined WARPX_DIM_1D_Z) + const ParticleReal p_x = 0.0; + const ParticleReal p_y = 0.0; +#elif (defined WARPX_DIM_RZ) const ParticleReal p_theta = p.rdata(PIdx::theta); const ParticleReal p_x = p_pos0*std::cos(p_theta); const ParticleReal p_y = p_pos0*std::sin(p_theta); +#elif (defined WARPX_DIM_XZ) + const ParticleReal p_x = p_pos0; + const ParticleReal p_y = 0.0; #else const ParticleReal p_pos1 = p.pos(1); const ParticleReal p_x = p_pos0; @@ -351,6 +386,20 @@ void BeamRelevant::ComputeDiags (int step) m_data[13] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c; m_data[14] = charge; amrex::ignore_unused(y_mean, y_ms, yuy); +#elif (defined WARPX_DIM_1D_Z) + m_data[0] = z_mean; + m_data[1] = ux_mean * m; + m_data[2] = uy_mean * m; + m_data[3] = uz_mean * m; + m_data[4] = gm_mean; + m_data[5] = std::sqrt(z_ms); + m_data[6] = std::sqrt(ux_ms) * m; + m_data[7] = std::sqrt(uy_ms) * m; + m_data[8] = std::sqrt(uz_ms) * m; + m_data[9] = std::sqrt(gm_ms); + m_data[10] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c; + m_data[11] = charge; + amrex::ignore_unused(x_mean, x_ms, xux, y_mean, y_ms, yuy); #endif } // end loop over species diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp index 378a11989..d83669ce5 100644 --- a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp @@ -99,7 +99,9 @@ void FieldEnergy::ComputeDiags (int step) // get cell size Geometry const & geom = warpx.Geom(lev); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + auto dV = geom.CellSize(0); +#elif (AMREX_SPACEDIM == 2) auto dV = geom.CellSize(0) * geom.CellSize(1); #elif (AMREX_SPACEDIM == 3) auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); diff --git a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp index c835feafd..1a0c9029d 100644 --- a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp @@ -186,7 +186,9 @@ void FieldMomentum::ComputeDiags (int step) // Get cell size amrex::Geometry const & geom = warpx.Geom(lev); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + auto dV = geom.CellSize(0); +#elif (AMREX_SPACEDIM == 2) auto dV = geom.CellSize(0) * geom.CellSize(1); #elif (AMREX_SPACEDIM == 3) auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index 648617719..adaa3d722 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -189,7 +189,9 @@ bool FieldProbe::ProbeInDomain () const * value in a position array will be the y value, but in the case of 2D, prob_lo[1] * and prob_hi[1] refer to z. This is a result of warpx.Geom(lev). */ -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + return z_probe >= prob_lo[1] && z_probe < prob_hi[1]; +#elif (AMREX_SPACEDIM == 2) return x_probe >= prob_lo[0] && x_probe < prob_hi[0] && z_probe >= prob_lo[1] && z_probe < prob_hi[1]; #else @@ -251,7 +253,10 @@ void FieldProbe::ComputeDiags (int step) { const auto cell_size = gm.CellSizeArray(); const int i_probe = static_cast<int>(amrex::Math::floor((x_probe - prob_lo[0]) / cell_size[0])); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const int j_probe = 0; + const int k_probe = 0; +#elif (AMREX_SPACEDIM == 2) const int j_probe = static_cast<int>(amrex::Math::floor((z_probe - prob_lo[1]) / cell_size[1])); const int k_probe = 0; #elif(AMREX_SPACEDIM == 3) diff --git a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp index c30d5820f..f092d53e3 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp @@ -185,6 +185,8 @@ void ParticleExtrema::ComputeDiags (int step) int const index_z = 2; #elif (defined WARPX_DIM_XZ || defined WARPX_DIM_RZ) int const index_z = 1; +#elif (defined WARPX_DIM_1D_Z) + int const index_z = 0; #endif // loop over species @@ -211,6 +213,8 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::cos(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMin(xmin); +#elif (defined WARPX_DIM_1D_Z) + Real xmin = 0.0_rt; #else Real xmin = ReduceMin( myspc, [=] AMREX_GPU_HOST_DEVICE (const PType& p) @@ -224,6 +228,8 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::cos(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMax(xmax); +#elif (defined WARPX_DIM_1D_Z) + Real xmax = 0.0_rt; #else Real xmax = ReduceMax( myspc, [=] AMREX_GPU_HOST_DEVICE (const PType& p) @@ -237,7 +243,7 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::sin(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMin(ymin); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) Real ymin = 0.0_rt; #else Real ymin = ReduceMin( myspc, @@ -252,7 +258,7 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::sin(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMax(ymax); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) Real ymax = 0.0_rt; #else Real ymax = ReduceMax( myspc, diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 392ef81b5..08d12dde6 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -38,7 +38,9 @@ WarpX::ComputeDt () if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + deltat = cfl * dx[0] / PhysConst::c; +#elif (AMREX_SPACEDIM == 2) deltat = cfl * std::min(dx[0], dx[1]) / PhysConst::c; #else deltat = cfl * std::min(dx[0], std::min(dx[1], dx[2])) / PhysConst::c; @@ -84,10 +86,13 @@ WarpX::PrintDtDxDyDz () for (int lev=0; lev <= max_level; lev++) { const amrex::Real* dx_lev = geom[lev].CellSize(); amrex::Print() << "Level " << lev << ": dt = " << dt[lev] +#if (defined WARPX_DIM_1D_Z) + << " ; dz = " << dx_lev[0] << '\n'; +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) << " ; dx = " << dx_lev[0] -#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) << " ; dz = " << dx_lev[1] << '\n'; #elif (defined WARPX_DIM_3D) + << " ; dx = " << dx_lev[0] << " ; dy = " << dx_lev[1] << " ; dz = " << dx_lev[2] << '\n'; #endif diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index d04a54ae3..770887cfe 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -253,7 +253,6 @@ WarpX::Evolve (int numsteps) bool move_j = is_synchronized; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. - int num_moved = MoveWindow(step+1, move_j); mypc->ContinuousFluxInjection(dt[0]); diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 5f33bbfb4..582890aa8 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -400,7 +400,9 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab> // Set the value of beta amrex::Array<amrex::Real,AMREX_SPACEDIM> beta_solver = -# if (AMREX_SPACEDIM==2) +# if (AMREX_SPACEDIM==1) + {{ beta[2] }}; // beta_x and beta_z +# elif (AMREX_SPACEDIM==2) {{ beta[0], beta[2] }}; // beta_x and beta_z # else {{ beta[0], beta[1], beta[2] }}; @@ -463,7 +465,13 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab> if (do_electrostatic == ElectrostaticSolverAlgo::LabFrame) { for (int lev = 0; lev <= max_level; ++lev) { -#if (AMREX_SPACEDIM==2) +#if (AMREX_SPACEDIM==1) + mlmg.getGradSolution( + {amrex::Array<amrex::MultiFab*,1>{ + get_pointer_Efield_fp(lev, 2) + }} + ); +#elif (AMREX_SPACEDIM==2) mlmg.getGradSolution( {amrex::Array<amrex::MultiFab*,2>{ get_pointer_Efield_fp(lev, 0),get_pointer_Efield_fp(lev, 2) @@ -578,21 +586,28 @@ WarpX::computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > #endif for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - const Real inv_dx = 1._rt/dx[0]; #if (AMREX_SPACEDIM == 3) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dy = 1._rt/dx[1]; const Real inv_dz = 1._rt/dx[2]; -#else +#elif (AMREX_SPACEDIM == 2) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dz = 1._rt/dx[1]; +#else + const Real inv_dz = 1._rt/dx[0]; #endif +#if (AMREX_SPACEDIM >= 2) const Box& tbx = mfi.tilebox( E[lev][0]->ixType().toIntVect() ); +#endif #if (AMREX_SPACEDIM == 3) const Box& tby = mfi.tilebox( E[lev][1]->ixType().toIntVect() ); #endif const Box& tbz = mfi.tilebox( E[lev][2]->ixType().toIntVect() ); const auto& phi_arr = phi[lev]->array(mfi); +#if (AMREX_SPACEDIM >= 2) const auto& Ex_arr = (*E[lev][0])[mfi].array(); +#endif #if (AMREX_SPACEDIM == 3) const auto& Ey_arr = (*E[lev][1])[mfi].array(); #endif @@ -631,7 +646,7 @@ WarpX::computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > +(beta_y*beta_z-1)*inv_dz*( phi_arr(i,j,k+1)-phi_arr(i,j,k) ); } ); -#else +#elif (AMREX_SPACEDIM == 2) amrex::ParallelFor( tbx, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Ex_arr(i,j,k) += @@ -646,6 +661,14 @@ WarpX::computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > +(beta_y*beta_z-1)*inv_dz*( phi_arr(i,j+1,k)-phi_arr(i,j,k) ); } ); +#else + amrex::ParallelFor( tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ez_arr(i,j,k) += + +(beta_y*beta_z-1)*inv_dz*( phi_arr(i+1,j,k)-phi_arr(i,j,k) ); + } + ); + ignore_unused(beta_x); #endif } } @@ -680,12 +703,15 @@ WarpX::computeB (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > #endif for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - const Real inv_dx = 1._rt/dx[0]; #if (AMREX_SPACEDIM == 3) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dy = 1._rt/dx[1]; const Real inv_dz = 1._rt/dx[2]; -#else +#elif (AMREX_SPACEDIM == 2) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dz = 1._rt/dx[1]; +#else + const Real inv_dz = 1._rt/dx[0]; #endif const Box& tbx = mfi.tilebox( B[lev][0]->ixType().toIntVect() ); const Box& tby = mfi.tilebox( B[lev][1]->ixType().toIntVect() ); @@ -728,7 +754,7 @@ WarpX::computeB (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > + phi_arr(i+1,j+1,k)-phi_arr(i,j+1,k))); } ); -#else +#elif (AMREX_SPACEDIM == 2) amrex::ParallelFor( tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Bx_arr(i,j,k) += inv_c * ( @@ -746,6 +772,18 @@ WarpX::computeB (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > +beta_y*inv_dx*( phi_arr(i+1,j,k)-phi_arr(i,j,k) )); } ); +#else + amrex::ParallelFor( tbx, tby, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Bx_arr(i,j,k) += inv_c * ( + -beta_y*inv_dz*( phi_arr(i+1,j,k)-phi_arr(i,j,k) )); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + By_arr(i,j,k) += inv_c * ( + +beta_x*inv_dz*(phi_arr(i+1,j,k)-phi_arr(i,j,k))); + } + ); + ignore_unused(beta_z,tbz,Bz_arr); #endif } } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp index 02b712679..46b32a8ff 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp @@ -165,9 +165,11 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( #else // Calculate relevant coefficients +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) amrex::Real const cdt_over_dx = PhysConst::c*dt*m_h_stencil_coefs_x[0]; amrex::Real const coef1_x = (1._rt - cdt_over_dx)/(1._rt + cdt_over_dx); amrex::Real const coef2_x = 2._rt*cdt_over_dx/(1._rt + cdt_over_dx) / PhysConst::c; +#endif #ifdef WARPX_DIM_3D amrex::Real const cdt_over_dy = PhysConst::c*dt*m_h_stencil_coefs_y[0]; amrex::Real const coef1_y = (1._rt - cdt_over_dy)/(1._rt + cdt_over_dy); @@ -177,8 +179,10 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( amrex::Real const coef1_z = (1._rt - cdt_over_dz)/(1._rt + cdt_over_dz); amrex::Real const coef2_z = 2._rt*cdt_over_dz/(1._rt + cdt_over_dz) / PhysConst::c; +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) bool const apply_lo_x = (field_boundary_lo[0] == FieldBoundaryType::Absorbing_SilverMueller); bool const apply_hi_x = (field_boundary_hi[0] == FieldBoundaryType::Absorbing_SilverMueller); +#endif #ifdef WARPX_DIM_3D bool const apply_lo_y = (field_boundary_lo[1] == FieldBoundaryType::Absorbing_SilverMueller); bool const apply_hi_y = (field_boundary_hi[1] == FieldBoundaryType::Absorbing_SilverMueller); @@ -201,10 +205,14 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // Extract field data for this grid/tile Array4<Real> const& Ex = Efield[0]->array(mfi); Array4<Real> const& Ey = Efield[1]->array(mfi); +#ifndef WARPX_DIM_1D_Z Array4<Real> const& Ez = Efield[2]->array(mfi); +#endif Array4<Real> const& Bx = Bfield[0]->array(mfi); Array4<Real> const& By = Bfield[1]->array(mfi); +#ifndef WARPX_DIM_1D_Z Array4<Real> const& Bz = Bfield[2]->array(mfi); +#endif // Extract the tileboxes for which to loop Box tbx = mfi.tilebox(Bfield[0]->ixType().toIntVect()); @@ -244,18 +252,27 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // At the -z boundary (innermost guard cell) if ( apply_lo_z && ( j==domain_box.smallEnd(1)-1 ) ) Bx(i,j,k) = coef1_z * Bx(i,j,k) + coef2_z * Ey(i,j+1,k); +#elif WARPX_DIM_1D_Z + // At the +z boundary (innermost guard cell) + if ( apply_hi_z && ( i==domain_box.bigEnd(0)+1 ) ) + Bx(i,j,k) = coef1_z * Bx(i,j,k) - coef2_z * Ey(i,j,k); + // At the -z boundary (innermost guard cell) + if ( apply_lo_z && ( i==domain_box.smallEnd(0)-1 ) ) + Bx(i,j,k) = coef1_z * Bx(i,j,k) + coef2_z * Ey(i+1,j,k); #endif }, // Apply Boundary condition to By [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) // At the +x boundary (innermost guard cell) if ( apply_hi_x && ( i==domain_box.bigEnd(0)+1 ) ) By(i,j,k) = coef1_x * By(i,j,k) - coef2_x * Ez(i,j,k); // At the -x boundary (innermost guard cell) if ( apply_lo_x && ( i==domain_box.smallEnd(0)-1 ) ) By(i,j,k) = coef1_x * By(i,j,k) + coef2_x * Ez(i+1,j,k); +#endif #ifdef WARPX_DIM_3D // At the +z boundary (innermost guard cell) if ( apply_hi_z && ( k==domain_box.bigEnd(2)+1 ) ) @@ -270,18 +287,27 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // At the -z boundary (innermost guard cell) if ( apply_lo_z && ( j==domain_box.smallEnd(1)-1 ) ) By(i,j,k) = coef1_z * By(i,j,k) - coef2_z * Ex(i,j+1,k); +#elif WARPX_DIM_1D_Z + // At the +z boundary (innermost guard cell) + if ( apply_hi_z && ( i==domain_box.bigEnd(0)+1 ) ) + By(i,j,k) = coef1_z * By(i,j,k) + coef2_z * Ex(i,j,k); + // At the -z boundary (innermost guard cell) + if ( apply_lo_z && ( i==domain_box.smallEnd(0)-1 ) ) + By(i,j,k) = coef1_z * By(i,j,k) - coef2_z * Ex(i+1,j,k); #endif }, // Apply Boundary condition to Bz [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) // At the +x boundary (innermost guard cell) if ( apply_hi_x && ( i==domain_box.bigEnd(0)+1 ) ) Bz(i,j,k) = coef1_x * Bz(i,j,k) + coef2_x * Ey(i,j,k); // At the -x boundary (innermost guard cell) if ( apply_lo_x && ( i==domain_box.smallEnd(0)-1 ) ) Bz(i,j,k) = coef1_x * Bz(i,j,k) - coef2_x * Ey(i+1,j,k); +#endif #ifdef WARPX_DIM_3D // At the +y boundary (innermost guard cell) if ( apply_hi_y && ( j==domain_box.bigEnd(1)+1 ) ) @@ -289,6 +315,8 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // At the -y boundary (innermost guard cell) if ( apply_lo_y && ( j==domain_box.smallEnd(1)-1 ) ) Bz(i,j,k) = coef1_y * Bz(i,j,k) + coef2_y * Ex(i,j+1,k); +#elif WARPX_DIM_1D_Z + ignore_unused(i,j,k); #endif } ); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index ac19984cd..318081159 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -176,7 +176,6 @@ void FiniteDifferenceSolver::EvolveECartesian ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k) - PhysConst::mu0 * jy(i, j, k) ); - }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H index 216723553..d7946c97e 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H @@ -71,6 +71,13 @@ struct CartesianCKCAlgorithm { constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt; constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt; constexpr Real alphay=0._rt; +#elif defined WARPX_DIM_1D_Z + Real const alphaz = inv_dz; + // Other coefficients are 0 in 1D Cartesian + // (and will actually not be used in the stencil) + constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt; + constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt, betaxz=0._rt, betazx=0._rt; + constexpr Real alphay=0._rt, alphax=0._rt; #endif // Store the coefficients in array `stencil_coefs`, in prescribed order @@ -124,14 +131,16 @@ struct CartesianCKCAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) amrex::Real const alphax = coefs_x[1]; -#if defined WARPX_DIM_3D - amrex::Real const betaxy = coefs_x[2]; -#endif amrex::Real const betaxz = coefs_x[3]; +#endif #if defined WARPX_DIM_3D + amrex::Real const betaxy = coefs_x[2]; amrex::Real const gammax = coefs_x[4]; #endif + #if defined WARPX_DIM_3D return alphax * (F(i+1,j ,k ,ncomp) - F(i, j, k ,ncomp)) + betaxy * (F(i+1,j+1,k ,ncomp) - F(i ,j+1,k ,ncomp) @@ -146,6 +155,9 @@ struct CartesianCKCAlgorithm { return alphax * (F(i+1,j ,k ,ncomp) - F(i, j, k ,ncomp)) + betaxz * (F(i+1,j+1,k ,ncomp) - F(i ,j+1,k ,ncomp) + F(i+1,j-1,k ,ncomp) - F(i ,j-1,k ,ncomp)); +#elif (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 #endif } @@ -158,8 +170,14 @@ struct CartesianCKCAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else amrex::Real const inv_dx = coefs_x[0]; return inv_dx*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); +#endif } /** @@ -186,10 +204,10 @@ struct CartesianCKCAlgorithm { + F(i-1,j+1,k+1,ncomp) - F(i-1,j ,k+1,ncomp) + F(i+1,j+1,k-1,ncomp) - F(i+1,j ,k-1,ncomp) + F(i-1,j+1,k-1,ncomp) - F(i-1,j ,k-1,ncomp)); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #else amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); @@ -210,7 +228,7 @@ struct CartesianCKCAlgorithm { amrex::Real const inv_dy = coefs_y[0]; amrex::ignore_unused(n_coefs_y); return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); return 0._rt; // 2D Cartesian: derivative along y is 0 @@ -233,7 +251,9 @@ struct CartesianCKCAlgorithm { amrex::ignore_unused(n_coefs_z); Real const alphaz = coefs_z[1]; +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) Real const betazx = coefs_z[2]; +#endif #if defined WARPX_DIM_3D Real const betazy = coefs_z[3]; Real const gammaz = coefs_z[4]; @@ -252,6 +272,8 @@ struct CartesianCKCAlgorithm { return alphaz * (F(i ,j+1,k ,ncomp) - F(i ,j ,k ,ncomp)) + betazx * (F(i+1,j+1,k ,ncomp) - F(i+1,j ,k ,ncomp) + F(i-1,j+1,k ,ncomp) - F(i-1,j ,k ,ncomp)); +#elif (defined WARPX_DIM_1D_Z) + return alphaz * (F(i+1 ,j,k ,ncomp) - F(i ,j ,k ,ncomp)); #endif } @@ -269,6 +291,8 @@ struct CartesianCKCAlgorithm { return inv_dz*( F(i,j,k,ncomp) - F(i,j,k-1,ncomp) ); #elif (defined WARPX_DIM_XZ) return inv_dz*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); #endif } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H index 5c1d68687..01d83804b 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H @@ -74,8 +74,13 @@ struct CartesianNodalAlgorithm { int const i, int const j, int const k, int const ncomp=0 ) { using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + ignore_unused(i, j, k, coefs_x, ncomp, F); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else Real const inv_dx = coefs_x[0]; return 0.5_rt*inv_dx*( F(i+1,j,k,ncomp) - F(i-1,j,k,ncomp) ); +#endif } /** @@ -88,8 +93,14 @@ struct CartesianNodalAlgorithm { amrex::Real const * const coefs_x, int const n_coefs_x, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + ignore_unused(i, j, k, coefs_x, n_coefs_x, ncomp, F); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else return UpwardDx( F, coefs_x, n_coefs_x, i, j, k ,ncomp); // For CartesianNodalAlgorithm, UpwardDx and DownwardDx are equivalent +#endif } /** @@ -106,9 +117,9 @@ struct CartesianNodalAlgorithm { #if defined WARPX_DIM_3D Real const inv_dy = coefs_y[0]; return 0.5_rt*inv_dy*( F(i,j+1,k,ncomp) - F(i,j-1,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_1D_Z) ignore_unused(i, j, k, coefs_y, ncomp, F); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #endif } @@ -142,6 +153,8 @@ struct CartesianNodalAlgorithm { return 0.5_rt*inv_dz*( F(i,j,k+1,ncomp) - F(i,j,k-1,ncomp) ); #elif (defined WARPX_DIM_XZ) return 0.5_rt*inv_dz*( F(i,j+1,k,ncomp) - F(i,j-1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return 0.5_rt*inv_dz*( F(i+1,j,k,ncomp) - F(i-1,j,k,ncomp) ); #endif } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H index 4fe813754..f2c8553d5 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -71,8 +71,14 @@ struct CartesianYeeAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else amrex::Real const inv_dx = coefs_x[0]; return inv_dx*( F(i+1,j,k,ncomp) - F(i,j,k,ncomp) ); +#endif } /** @@ -84,8 +90,14 @@ struct CartesianYeeAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else amrex::Real const inv_dx = coefs_x[0]; return inv_dx*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); +#endif } /** @@ -101,10 +113,10 @@ struct CartesianYeeAlgorithm { Real const inv_dy = coefs_y[0]; amrex::ignore_unused(n_coefs_y); return inv_dy*( F(i,j+1,k,ncomp) - F(i,j,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #else amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); @@ -125,10 +137,10 @@ struct CartesianYeeAlgorithm { Real const inv_dy = coefs_y[0]; amrex::ignore_unused(n_coefs_y); return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #else amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); @@ -149,6 +161,8 @@ struct CartesianYeeAlgorithm { return inv_dz*( F(i,j,k+1,ncomp) - F(i,j,k,ncomp) ); #elif (defined WARPX_DIM_XZ) return inv_dz*( F(i,j+1,k,ncomp) - F(i,j,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz*( F(i+1,j,k,ncomp) - F(i,j,k,ncomp) ); #endif } @@ -167,6 +181,8 @@ struct CartesianYeeAlgorithm { return inv_dz*( F(i,j,k,ncomp) - F(i,j,k-1,ncomp) ); #elif (defined WARPX_DIM_XZ) return inv_dz*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); #endif } diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 2ac27227d..6f1ebad78 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -32,6 +32,10 @@ void warpx_computedivb(int i, int j, int k, int dcomp, divB(i,j,0,dcomp) = (Bx(i+1,j ,0) - Bx(i,j,0))*dxinv + (Bz(i ,j+1,0) - Bz(i,j,0))*dzinv; amrex::ignore_unused(k, By, dyinv); +#elif defined WARPX_DIM_1D_Z + divB(i,0,0,dcomp) = (Bz(i+1,0 ,0) - Bz(i,0,0))*dzinv; + amrex::ignore_unused(j, Bx, dxinv); + amrex::ignore_unused(k, By, dyinv); #elif defined WARPX_DIM_RZ const amrex::Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); const amrex::Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index 8a16e2a9b..601bd9ffc 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -80,9 +80,17 @@ void BilinearFilter::ComputeStencils(){ stencil_z.resize( 1u + npass_each_dir[1] ); compute_stencil(stencil_x, npass_each_dir[0]); compute_stencil(stencil_z, npass_each_dir[1]); +#elif (AMREX_SPACEDIM == 1) + // npass_each_dir = npass_z + stencil_z.resize( 1u + npass_each_dir[0] ); + compute_stencil(stencil_z, npass_each_dir[0]); #endif slen = stencil_length_each_dir.dim3(); #if (AMREX_SPACEDIM == 2) slen.z = 1; #endif +#if (AMREX_SPACEDIM == 1) + slen.y = 1; + slen.z = 1; +#endif } diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 654a1aa56..b72d6cec6 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -92,7 +92,9 @@ void Filter::DoFilter (const Box& tbx, Array4<Real > const& dst, int scomp, int dcomp, int ncomp) { +#if (AMREX_SPACEDIM >= 2) amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); +#endif #if (AMREX_SPACEDIM == 3) amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); #endif @@ -129,7 +131,7 @@ void Filter::DoFilter (const Box& tbx, dst(i,j,k,dcomp+n) = d; }); -#else +#elif (AMREX_SPACEDIM == 2) AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, { Real d = 0.0; @@ -155,6 +157,32 @@ void Filter::DoFilter (const Box& tbx, dst(i,j,k,dcomp+n) = d; }); +#elif (AMREX_SPACEDIM == 1) + AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, + { + Real d = 0.0; + + // Pad source array with zeros beyond ghost cells + // for out-of-bound accesses due to large-stencil operations + const auto src_zeropad = [src] (const int jj, const int kk, const int ll, const int nn) noexcept + { + return src.contains(jj,kk,ll) ? src(jj,kk,ll,nn) : 0.0_rt; + }; + + for (int iz=0; iz < slen_local.z; ++iz){ + for (int iy=0; iy < slen_local.y; ++iy){ + for (int ix=0; ix < slen_local.x; ++ix){ + Real sss = sz[iy]; + d += sss*( src_zeropad(i-ix,j,k,scomp+n) + +src_zeropad(i+ix,j,k,scomp+n)); + } + } + } + + dst(i,j,k,dcomp+n) = d; + }); +#else + amrex::Abort("Filter not implemented for the current geometry!"); #endif } @@ -247,7 +275,9 @@ void Filter::DoFilter (const Box& tbx, const auto lo = amrex::lbound(tbx); const auto hi = amrex::ubound(tbx); // tmp and dst are of type Array4 (Fortran ordering) +#if (AMREX_SPACEDIM >= 2) amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); +#endif #if (AMREX_SPACEDIM == 3) amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); #endif @@ -267,8 +297,10 @@ void Filter::DoFilter (const Box& tbx, for (int ix=0; ix < slen.x; ++ix){ #if (AMREX_SPACEDIM == 3) Real sss = sx[ix]*sy[iy]*sz[iz]; -#else +#elif (AMREX_SPACEDIM == 2) Real sss = sx[ix]*sz[iy]; +#else + Real sss = sz[ix]; #endif // 3 nested loop on 3D array for (int k = lo.z; k <= hi.z; ++k) { @@ -284,11 +316,16 @@ void Filter::DoFilter (const Box& tbx, +tmp(i+ix,j-iy,k+iz,scomp+n) +tmp(i-ix,j+iy,k+iz,scomp+n) +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else +#elif (AMREX_SPACEDIM == 2) dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) +tmp(i+ix,j-iy,k,scomp+n) +tmp(i-ix,j+iy,k,scomp+n) +tmp(i+ix,j+iy,k,scomp+n)); +#elif (AMREX_SPACEDIM == 1) + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j,k,scomp+n) + +tmp(i+ix,j,k,scomp+n)); +#else + amrex::Abort("Filter not implemented for the current geometry!"); #endif } } diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 224bfb946..fa88eb854 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -23,6 +23,7 @@ using namespace amrex; NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set, amrex::Real cdtodz, bool nodal_gather){ +#if (AMREX_SPACEDIM >= 2) // Store parameters into class data members m_coeff_set = coeff_set; m_cdtodz = cdtodz; @@ -35,10 +36,15 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set, amrex::Real cdto stencil_length_each_dir = {1,5}; slen = {1,5,1}; #endif +#else + amrex::ignore_unused(coeff_set, cdtodz, nodal_gather); + amrex::Abort("NCIGodfreyFilter not implemented in 1D!"); +#endif } void NCIGodfreyFilter::ComputeStencils(){ +#if (AMREX_SPACEDIM >= 2) using namespace warpx::nci_godfrey; // Sanity checks: filter length shoulz be 5 in z @@ -128,4 +134,7 @@ void NCIGodfreyFilter::ComputeStencils(){ Gpu::copyAsync(Gpu::hostToDevice,h_stencil_z.begin(),h_stencil_z.end(),stencil_z.begin()); Gpu::synchronize(); +#else + amrex::Abort("NCIGodfreyFilter not implemented in 1D!"); +#endif } diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index af6a7a208..5710cab9a 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -81,10 +81,16 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) // NOTE: When periodic boundaries are used, default injection range is set to mother grid dimensions. const amrex::Geometry& geom = WarpX::GetInstance().Geom(0); if( geom.isPeriodic(0) ) { +# ifndef WARPX_DIM_1D_Z xmin = geom.ProbLo(0); xmax = geom.ProbHi(0); +# else + zmin = geom.ProbLo(0); + zmax = geom.ProbHi(0); +# endif } +# ifndef WARPX_DIM_1D_Z if( geom.isPeriodic(1) ) { # ifndef WARPX_DIM_3D zmin = geom.ProbLo(1); @@ -94,6 +100,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) ymax = geom.ProbHi(1); # endif } +# endif # ifdef WARPX_DIM_3D if( geom.isPeriodic(2) ) { @@ -246,16 +253,18 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) flux_normal_axis = 0; } #else +# ifndef WARPX_DIM_1D_Z if (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") { flux_normal_axis = 0; } +# endif #endif #ifdef WARPX_DIM_3D - else if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") { + if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") { flux_normal_axis = 1; } #endif - else if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") { + if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") { flux_normal_axis = AMREX_SPACEDIM-1; } #ifdef WARPX_DIM_3D @@ -263,8 +272,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) #else # ifdef WARPX_DIM_RZ std::string flux_normal_axis_help = "'r' or 'z'."; -# else +# elif WARPX_DIM_XZ std::string flux_normal_axis_help = "'x' or 'z'."; +# else + std::string flux_normal_axis_help = "'z'."; # endif #endif AMREX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0, @@ -282,7 +293,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) } else if (injection_style == "nuniformpercell") { // Note that for RZ, three numbers are expected, r, theta, and z. // For 2D, only two are expected. The third is overwritten with 1. -#if WARPX_DIM_XZ + // For 1D, only one is expected. The second and third are overwritten with 1. +#if defined(WARPX_DIM_1D_Z) + constexpr int num_required_ppc_each_dim = 1; +#elif defined(WARPX_DIM_XZ) constexpr int num_required_ppc_each_dim = 2; #else constexpr int num_required_ppc_each_dim = 3; @@ -292,6 +306,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) #if WARPX_DIM_XZ num_particles_per_cell_each_dim.push_back(1); #endif +#if WARPX_DIM_1D_Z + num_particles_per_cell_each_dim.push_back(1); // overwrite 2nd number with 1 + num_particles_per_cell_each_dim.push_back(1); // overwrite 3rd number with 1 +#endif #if WARPX_DIM_RZ if (WarpX::n_rz_azimuthal_modes > 1) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 41fda5567..f92fe80aa 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -86,11 +86,13 @@ WarpX::PostProcessBaseGrids (BoxArray& ba0) const klo += domlo[2]; khi += domlo[2]; #endif +#if (AMREX_SPACEDIM >= 2) for (int j = 0; j < numprocs[1]; ++j) { int jlo = (j < extra[1]) ? j*(sz[1]+1) : (j*sz[1]+extra[1]); int jhi = (j < extra[1]) ? jlo+(sz[1]+1)-1 : jlo+sz[1]-1; jlo += domlo[1]; jhi += domlo[1]; +#endif for (int i = 0; i < numprocs[0]; ++i) { int ilo = (i < extra[0]) ? i*(sz[0]+1) : (i*sz[0]+extra[0]); int ihi = (i < extra[0]) ? ilo+(sz[0]+1)-1 : ilo+sz[0]-1; @@ -358,6 +360,7 @@ WarpX::computeMaxStepBoostAccelerator(const amrex::Geometry& a_geom){ void WarpX::InitNCICorrector () { +#if !(defined WARPX_DIM_1D_Z) if (WarpX::use_fdtd_nci_corr) { for (int lev = 0; lev <= max_level; ++lev) @@ -367,8 +370,10 @@ WarpX::InitNCICorrector () amrex::Real dz, cdtodz; if (AMREX_SPACEDIM == 3){ dz = dx[2]; - }else{ + }else if(AMREX_SPACEDIM == 2){ dz = dx[1]; + }else{ + dz = dx[0]; } cdtodz = PhysConst::c * dt[lev] / dz; @@ -385,6 +390,7 @@ WarpX::InitNCICorrector () nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); } } +#endif } void @@ -663,13 +669,20 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( #endif // Shift required in the x-, y-, or z- position // depending on the index type of the multifab +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0._rt; amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; #else + amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1._rt - x_nodal_flag[2]) * dx_lev[2] * 0.5_rt; @@ -682,13 +695,20 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( #ifdef AMREX_USE_EB if(geom_data_y(i, j, k)<=0) return; #endif +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0._rt; amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; #elif (AMREX_SPACEDIM==3) + amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1._rt - y_nodal_flag[2]) * dx_lev[2] * 0.5_rt; @@ -701,13 +721,20 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( #ifdef AMREX_USE_EB if(geom_data_z(i, j, k)<=0) return; #endif +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0._rt; amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; #elif (AMREX_SPACEDIM==3) + amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1._rt - z_nodal_flag[2]) * dx_lev[2] * 0.5_rt; diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 5df0635e7..e974e5764 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -58,7 +58,9 @@ include $(AMREX_HOME)/Tools/GNUMake/Make.defs ifeq ($(DIM),3) DEFINES += -DWARPX_DIM_3D -else +endif + +ifeq ($(DIM),2) ifeq ($(USE_RZ),TRUE) DEFINES += -DWARPX_DIM_RZ else @@ -66,6 +68,10 @@ else endif endif +ifeq ($(DIM),1) + DEFINES += -DWARPX_DIM_1D_Z +endif + -include Make.package include $(WARPX_HOME)/Source/Make.package include $(WARPX_HOME)/Source/BoundaryConditions/Make.package diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp index 44848c915..b72e13cdc 100644 --- a/Source/Parallelization/GuardCellManager.cpp +++ b/Source/Parallelization/GuardCellManager.cpp @@ -115,6 +115,9 @@ guardCellManager::Init ( #elif (AMREX_SPACEDIM == 2) ng_alloc_EB = IntVect(ngx,ngz); ng_alloc_J = IntVect(ngJx,ngJz); +#elif (AMREX_SPACEDIM == 1) + ng_alloc_EB = IntVect(ngz); + ng_alloc_J = IntVect(ngJz); #endif // TODO Adding one cell for rho should not be necessary, given that the number of guard cells @@ -195,6 +198,8 @@ guardCellManager::Init ( IntVect ngFFT = IntVect(ngFFt_x, ngFFt_y, ngFFt_z); #elif (AMREX_SPACEDIM == 2) IntVect ngFFT = IntVect(ngFFt_x, ngFFt_z); +#elif (AMREX_SPACEDIM == 1) + IntVect ngFFT = IntVect(ngFFt_z); #endif // All boxes should have the same number of guard cells, to avoid temporary parallel copies: diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 2d1e1f04a..ecaddb570 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -979,7 +979,9 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) IntVect ng_depos_J = get_ng_depos_J(); if (WarpX::do_current_centering) { -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + ng_depos_J[0] += WarpX::current_centering_noz / 2; +#elif (AMREX_SPACEDIM == 2) ng_depos_J[0] += WarpX::current_centering_nox / 2; ng_depos_J[1] += WarpX::current_centering_noz / 2; #elif (AMREX_SPACEDIM == 3) diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index 04226e971..5995243f5 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -20,17 +20,17 @@ void warpx_interp (int j, int k, int l, { using namespace amrex; - // NOTE Indices (j,k,l) in the following refer to (x,z,-) in 2D and (x,y,z) in 3D + // NOTE Indices (j,k,l) in the following refer to (z,-,-) in 1D, (x,z,-) in 2D, and (x,y,z) in 3D // Refinement ratio const int rj = rr[0]; - const int rk = rr[1]; - const int rl = (AMREX_SPACEDIM == 2) ? 1 : rr[2]; + const int rk = (AMREX_SPACEDIM == 1) ? 1 : rr[1]; + const int rl = (AMREX_SPACEDIM <= 2) ? 1 : rr[2]; // Staggering (0: cell-centered; 1: nodal) const int sj = arr_stag[0]; - const int sk = arr_stag[1]; - const int sl = (AMREX_SPACEDIM == 2) ? 0 : arr_stag[2]; + const int sk = (AMREX_SPACEDIM == 1) ? 0 : arr_stag[1]; + const int sl = (AMREX_SPACEDIM <= 2) ? 0 : arr_stag[2]; // Number of points used for interpolation from coarse grid to fine grid const int nj = (sj == 0) ? 1 : 2; @@ -87,7 +87,21 @@ void warpx_interp_nd_bfield_x (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real bg = owx * Bxg(jg ,0,0) + + wx * Bxg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real bc = owx * Bxc(jg ,0,0) + + wx * Bxc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5_rt*(Bxf_zeropad(j-1,0,0) + Bxf_zeropad(j,0,0)); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real bg = owx * owy * Bxg(jg ,kg ,0) @@ -163,7 +177,21 @@ void warpx_interp_nd_bfield_y (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real bg = owx * Byg(jg ,0,0) + + wx * Byg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real bc = owx * Byc(jg ,0,0) + + wx * Byc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5_rt*(Byf_zeropad(j-1,0,0) + Byf_zeropad(j,0,0)); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real bg = owx * owy * Byg(jg ,kg ,0) @@ -241,7 +269,21 @@ void warpx_interp_nd_bfield_z (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real bg = owx * Bzg(jg ,0,0) + + wx * Bzg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real bc = owx * Bzc(jg ,0,0) + + wx * Bzc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real bf = Bzf_zeropad(j,0,0); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real bg = owx * owy * Bzg(jg ,kg ,0) @@ -318,7 +360,21 @@ void warpx_interp_nd_efield_x (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real eg = owx * Exg(jg ,0,0) + + wx * Exg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real ec = owx * Exc(jg ,0,0) + + wx * Exc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real ef = Exf_zeropad(j,0,0); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real eg = owx * owy * Exg(jg ,kg ,0) @@ -394,7 +450,23 @@ void warpx_interp_nd_efield_y (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) + + +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real eg = owx * Eyg(jg ,0,0) + + wx * Eyg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real ec = owx * Eyc(jg ,0,0) + + wx * Eyc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real ef = Eyf_zeropad(j,0,0); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real eg = owx * owy * Eyg(jg ,kg ,0) @@ -469,7 +541,21 @@ void warpx_interp_nd_efield_z (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real eg = owx * Ezg(jg ,0,0) + + wx * Ezg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real ec = owx * Ezc(jg ,0,0) + + wx * Ezc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real ef = 0.5_rt*(Ezf_zeropad(j-1,0,0) + Ezf_zeropad(j,0,0)); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real eg = owx * owy * Ezg(jg ,kg ,0) @@ -568,6 +654,9 @@ void warpx_interp (const int j, }; // Avoid compiler warnings +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(nox, noy, stencil_coeffs_x, stencil_coeffs_y); +#endif #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(noy, stencil_coeffs_y); #endif @@ -586,32 +675,43 @@ void warpx_interp (const int j, // Staggering (s = 0 if cell-centered, s = 1 if nodal) const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0]; +#if (AMREX_SPACEDIM >= 2) const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1]; +#endif #if (AMREX_SPACEDIM == 3) const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2]; #endif // Interpolate along j,k,l only if source MultiFab is staggered along j,k,l const bool interp_j = (sj == 0); +#if (AMREX_SPACEDIM >= 2) const bool interp_k = (sk == 0); +#endif #if (AMREX_SPACEDIM == 3) const bool interp_l = (sl == 0); #endif +#if (AMREX_SPACEDIM == 1) + const int noj = noz; +#elif (AMREX_SPACEDIM == 2) const int noj = nox; -#if (AMREX_SPACEDIM == 2) const int nok = noz; #elif (AMREX_SPACEDIM == 3) + const int noj = nox; const int nok = noy; const int nol = noz; #endif // Additional normalization factor const amrex::Real wj = (interp_j) ? 0.5_rt : 1.0_rt; +#if (AMREX_SPACEDIM == 1) + constexpr amrex::Real wk = 1.0_rt; + constexpr amrex::Real wl = 1.0_rt; +#elif (AMREX_SPACEDIM == 2) const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt; -#if (AMREX_SPACEDIM == 2) constexpr amrex::Real wl = 1.0_rt; #elif (AMREX_SPACEDIM == 3) + const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt; const amrex::Real wl = (interp_l) ? 0.5_rt : 1.0_rt; #endif @@ -620,11 +720,17 @@ void warpx_interp (const int j, const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j; // Min and max for interpolation loop along k +#if (AMREX_SPACEDIM == 1) + // k = 0 always + const int kmin = k; + const int kmax = k; +#else const int kmin = (interp_k) ? k - nok/2 + shift : k; const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k; +#endif // Min and max for interpolation loop along l -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM <= 2) // l = 0 always const int lmin = l; const int lmax = l; @@ -696,6 +802,11 @@ void warpx_interp (const int j, else // PSATD (finite-order interpolation) { +#if (AMREX_SPACEDIM == 1) + amrex::Abort("PSATD not implemented in 1D"); +#endif + +#if (AMREX_SPACEDIM >= 2) // 1D not implemented for PSATD amrex::Real cj = 1.0_rt; amrex::Real ck = 1.0_rt; amrex::Real cl = 1.0_rt; @@ -725,6 +836,7 @@ void warpx_interp (const int j, } } } +#endif //1D } dst_arr(j,k,l) = wj * wk * wl * res; diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index fa3a0effc..98da517d6 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -257,7 +257,9 @@ public: const amrex::Real dt = WarpX::GetInstance().getdt(lev); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); -#if defined WARPX_DIM_XZ +#if defined WARPX_DIM_1D_Z + auto dV = geom.CellSize(0); +#elif defined WARPX_DIM_XZ auto dV = geom.CellSize(0) * geom.CellSize(1); #elif defined WARPX_DIM_RZ amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box @@ -406,7 +408,9 @@ public: const amrex::Real dt = WarpX::GetInstance().getdt(lev); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); -#if defined WARPX_DIM_XZ +#if defined WARPX_DIM_1D_Z + auto dV = geom.CellSize(0); +#elif defined WARPX_DIM_XZ auto dV = geom.CellSize(0) * geom.CellSize(1); #elif defined WARPX_DIM_RZ amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index d651b1b6c..84db47ba5 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -58,16 +58,22 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) const bool do_ionization = ion_lev; - const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dzi = 1.0_rt/dx[2]; +#if (AMREX_SPACEDIM == 1) + const amrex::Real invvol = dzi; +#endif #if (AMREX_SPACEDIM == 2) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real invvol = dxi*dzi; #elif (defined WARPX_DIM_3D) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dyi = 1.0_rt/dx[1]; const amrex::Real invvol = dxi*dyi*dzi; #endif +#if (AMREX_SPACEDIM >= 2) const amrex::Real xmin = xyzmin[0]; +#endif #if (defined WARPX_DIM_3D) const amrex::Real ymin = xyzmin[1]; #endif @@ -105,6 +111,8 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, GetPosition(ip, xp, yp, zp); // --- Compute shape factors + Compute_shape_factor< depos_order > const compute_shape_factor; +#if (AMREX_SPACEDIM >= 2) // x direction // Get particle position in grid coordinates #if (defined WARPX_DIM_RZ) @@ -128,13 +136,12 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, // i: leftmost grid point that the particle touches amrex::Real sx[depos_order + 1] = {0._rt}; int i = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; if (rho_type[0] == NODE) { i = compute_shape_factor(sx, x); } else if (rho_type[0] == CELL) { i = compute_shape_factor(sx, x - 0.5_rt); } - +#endif //AMREX_SPACEDIM >= 2 #if (defined WARPX_DIM_3D) // y direction const amrex::Real y = (yp - ymin)*dyi; @@ -157,6 +164,13 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, } // Deposit charge into rho_arr +#if (defined WARPX_DIM_1D_Z) + for (int iz=0; iz<=depos_order; iz++){ + amrex::Gpu::Atomic::AddNoRet( + &rho_arr(lo.x+k+iz, 0, 0, 0), + sz[iz]*wq); + } +#endif #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 6c330e388..640ba4514 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -79,16 +79,22 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) const bool do_ionization = ion_lev; - const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dzi = 1.0_rt/dx[2]; +#if (AMREX_SPACEDIM == 1) + const amrex::Real invvol = dzi; +#endif #if (AMREX_SPACEDIM == 2) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real invvol = dxi*dzi; #elif (defined WARPX_DIM_3D) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dyi = 1.0_rt/dx[1]; const amrex::Real invvol = dxi*dyi*dzi; #endif +#if (AMREX_SPACEDIM >= 2) const amrex::Real xmin = xyzmin[0]; +#endif #if (defined WARPX_DIM_3D) const amrex::Real ymin = xyzmin[1]; #endif @@ -164,6 +170,8 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, const amrex::Real wqz = wq*invvol*vz; // --- Compute shape factors + Compute_shape_factor< depos_order > const compute_shape_factor; +#if (AMREX_SPACEDIM >= 2) // x direction // Get particle position after 1/2 push back in position #if (defined WARPX_DIM_RZ) @@ -181,7 +189,6 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, double sx_cell[depos_order + 1] = {0.}; int j_node = 0; int j_cell = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) { j_node = compute_shape_factor(sx_node, xmid); } @@ -202,6 +209,7 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell); int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell); int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell); +#endif //AMREX_SPACEDIM >= 2 #if (defined WARPX_DIM_3D) // y direction @@ -258,6 +266,19 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell); // Deposit current into jx_arr, jy_arr and jz_arr +#if (defined WARPX_DIM_1D_Z) + for (int iz=0; iz<=depos_order; iz++){ + amrex::Gpu::Atomic::AddNoRet( + &jx_arr(lo.x+l_jx+iz, 0, 0, 0), + sz_jx[iz]*wqx); + amrex::Gpu::Atomic::AddNoRet( + &jy_arr(lo.x+l_jy+iz, 0, 0, 0), + sz_jy[iz]*wqy); + amrex::Gpu::Atomic::AddNoRet( + &jz_arr(lo.x+l_jz+iz, 0, 0, 0), + sz_jz[iz]*wqz); + } +#endif #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ @@ -367,11 +388,15 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) bool const do_ionization = ion_lev; +#if !(defined WARPX_DIM_1D_Z) Real const dxi = 1.0_rt / dx[0]; -#if !(defined WARPX_DIM_RZ) +#endif +#if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) Real const dtsdx0 = dt*dxi; #endif +#if !(defined WARPX_DIM_1D_Z) Real const xmin = xyzmin[0]; +#endif #if (defined WARPX_DIM_3D) Real const dyi = 1.0_rt / dx[1]; Real const dtsdy0 = dt*dyi; @@ -389,6 +414,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const invdtdx = 1.0_rt / (dt*dx[2]); Real const invdtdz = 1.0_rt / (dt*dx[0]); Real const invvol = 1.0_rt / (dx[0]*dx[2]); +#elif (defined WARPX_DIM_1D_Z) + Real const invdtdz = 1.0_rt / (dt*dx[0]); + Real const invvol = 1.0_rt / (dx[2]); #endif #if (defined WARPX_DIM_RZ) @@ -396,8 +424,10 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, #endif Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); +#if !(defined WARPX_DIM_1D_Z) Real constexpr one_third = 1.0_rt / 3.0_rt; Real constexpr one_sixth = 1.0_rt / 6.0_rt; +#endif // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr #if defined(WARPX_USE_GPUCLOCK) @@ -429,7 +459,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, ParticleReal xp, yp, zp; GetPosition(ip, xp, yp, zp); +#if !(defined WARPX_DIM_1D_Z) Real const wqx = wq*invdtdx; +#endif #if (defined WARPX_DIM_3D) Real const wqy = wq*invdtdy; #endif @@ -476,9 +508,11 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, double const x_old = (rp_old - xmin)*dxi; #else // Keep these double to avoid bug in single precision +#if !(defined WARPX_DIM_1D_Z) double const x_new = (xp - xmin)*dxi; double const x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif +#endif #if (defined WARPX_DIM_3D) // Keep these double to avoid bug in single precision double const y_new = (yp - ymin)*dyi; @@ -492,6 +526,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; #elif (defined WARPX_DIM_XZ) Real const vy = uyp[ip]*gaminv; +#elif (defined WARPX_DIM_1D_Z) + Real const vx = uxp[ip]*gaminv; + Real const vy = uyp[ip]*gaminv; #endif // Shape factor arrays @@ -499,8 +536,10 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, // to possibly hold the factor for the old particle // which can be at a different grid location. // Keep these double to avoid bug in single precision +#if !(defined WARPX_DIM_1D_Z) double sx_new[depos_order + 3] = {0.}; double sx_old[depos_order + 3] = {0.}; +#endif #if (defined WARPX_DIM_3D) // Keep these double to avoid bug in single precision double sy_new[depos_order + 3] = {0.}; @@ -516,8 +555,10 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Compute_shape_factor< depos_order > compute_shape_factor; Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; +#if !(defined WARPX_DIM_1D_Z) const int i_new = compute_shape_factor(sx_new+1, x_new); const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new); +#endif #if (defined WARPX_DIM_3D) const int j_new = compute_shape_factor(sy_new+1, y_new); const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new); @@ -526,9 +567,11 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new); // computes min/max positions of current contributions +#if !(defined WARPX_DIM_1D_Z) int dil = 1, diu = 1; if (i_old < i_new) dil = 0; if (i_old > i_new) diu = 0; +#endif #if (defined WARPX_DIM_3D) int djl = 1, dju = 1; if (j_old < j_new) djl = 0; @@ -636,6 +679,23 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, #endif } } +#elif (defined WARPX_DIM_1D_Z) + + for (int k=dkl; k<=depos_order+2-dku; k++) { + amrex::Real sdxi = 0._rt; + sdxi += wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]); + amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+k_new-1+k, 0, 0, 0), sdxi); + } + for (int k=dkl; k<=depos_order+2-dku; k++) { + amrex::Real sdyj = 0._rt; + sdyj += wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]); + amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+k_new-1+k, 0, 0, 0), sdyj); + } + for (int k=dkl; k<=depos_order+1-dku; k++) { + amrex::Real sdzk = 0._rt; + sdzk += wqz*(sz_old[k] - sz_new[k]); + amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk); + } #endif } ); @@ -700,11 +760,18 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, amrex::Abort("Vay deposition not implemented in RZ geometry"); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(GetPosition, + wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab, + np_to_depose, dt, dx, xyzmin, lo, q, n_rz_azimuthal_modes); + amrex::Abort("Vay deposition not implemented in cartesian 1D geometry"); +#endif + #if !defined(AMREX_USE_GPU) amrex::ignore_unused(cost, load_balance_costs_update_algo); #endif -#if !(defined WARPX_DIM_RZ) +#if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) amrex::ignore_unused(n_rz_azimuthal_modes); // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1 @@ -1001,6 +1068,6 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, amrex::The_Managed_Arena()->free(cost_real); } # endif -#endif // #if !(defined WARPX_DIM_RZ) +#endif // #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) } #endif // CURRENTDEPOSITION_H_ diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 562f5b242..23a27f330 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -67,17 +67,25 @@ void doGatherShapeN (const amrex::ParticleReal xp, amrex::ignore_unused(yp); #endif +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(xp,yp); +#endif + #ifndef WARPX_DIM_RZ amrex::ignore_unused(n_rz_azimuthal_modes); #endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 3) const amrex::Real dxi = 1.0_rt/dx[0]; +#endif const amrex::Real dzi = 1.0_rt/dx[2]; #if (AMREX_SPACEDIM == 3) const amrex::Real dyi = 1.0_rt/dx[1]; #endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 3) const amrex::Real xmin = xyzmin[0]; +#endif #if (AMREX_SPACEDIM == 3) const amrex::Real ymin = xyzmin[1]; #endif @@ -88,6 +96,11 @@ void doGatherShapeN (const amrex::ParticleReal xp, constexpr int CELL = amrex::IndexType::CELL; // --- Compute shape factors + + Compute_shape_factor< depos_order > const compute_shape_factor; + Compute_shape_factor<depos_order - galerkin_interpolation > const compute_shape_factor_galerkin; + +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 3) // x direction // Get particle position #ifdef WARPX_DIM_RZ @@ -110,8 +123,6 @@ void doGatherShapeN (const amrex::ParticleReal xp, int j_cell = 0; int j_node_v = 0; int j_cell_v = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; - Compute_shape_factor<depos_order - galerkin_interpolation > const compute_shape_factor_galerkin; if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) { j_node = compute_shape_factor(sx_node, x); } @@ -136,6 +147,7 @@ void doGatherShapeN (const amrex::ParticleReal xp, int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell ); int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v); int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v); +#endif #if (AMREX_SPACEDIM == 3) // y direction @@ -214,7 +226,33 @@ void doGatherShapeN (const amrex::ParticleReal xp, // AMREX_SPACEDIM nested loops because the deposition // order can differ for each component of each field // when galerkin_interpolation is set to 1 -#if (AMREX_SPACEDIM == 2) + +#if (AMREX_SPACEDIM == 1) + // Gather field on particle Eyp from field on grid ey_arr + // Gather field on particle Exp from field on grid ex_arr + // Gather field on particle Bzp from field on grid bz_arr + for (int iz=0; iz<=depos_order; iz++){ + Eyp += sz_ey[iz]* + ey_arr(lo.x+l_ey+iz, 0, 0, 0); + Exp += sz_ex[iz]* + ex_arr(lo.x+l_ex+iz, 0, 0, 0); + Bzp += sz_bz[iz]* + bz_arr(lo.x+l_bz+iz, 0, 0, 0); + } + + // Gather field on particle Byp from field on grid by_arr + // Gather field on particle Ezp from field on grid ez_arr + // Gather field on particle Bxp from field on grid bx_arr + for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){ + Ezp += sz_ez[iz]* + ez_arr(lo.x+l_ez+iz, 0, 0, 0); + Bxp += sz_bx[iz]* + bx_arr(lo.x+l_bx+iz, 0, 0, 0); + Byp += sz_by[iz]* + by_arr(lo.x+l_by+iz, 0, 0, 0); + } + +#elif (AMREX_SPACEDIM == 2) // Gather field on particle Eyp from field on grid ey_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 559489f34..e7c47a6f4 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -134,6 +134,12 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec[1] == amrex::Real(0), "Laser propagation direction must be 0 along y in 2D"); #endif +#ifdef WARPX_DIM_1D_Z + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec[0] == amrex::Real(0), + "Laser propagation direction must be 0 along x in 1D"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec[1] == amrex::Real(0), + "Laser propagation direction must be 0 along y in 2D"); +#endif // Plane normal Real s = 1.0_rt / std::sqrt(m_nvec[0]*m_nvec[0] + m_nvec[1]*m_nvec[1] + m_nvec[2]*m_nvec[2]); @@ -168,9 +174,12 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) m_u_X = m_p_X; m_u_Y = m_p_Y; -#else +#elif (defined WARPX_DIM_XZ) m_u_X = CrossProduct({0., 1., 0.}, m_nvec); m_u_Y = {0., 1., 0.}; +#elif (defined WARPX_DIM_1D_Z) + m_u_X = {1., 0., 0.}; + m_u_Y = {0., 1., 0.}; #endif m_laser_injection_box= Geom(0).ProbDomain(); @@ -192,7 +201,10 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, // Sanity checks int dir = WarpX::moving_window_dir; std::vector<Real> windir(3, 0.0); -#if (AMREX_SPACEDIM==2) +#if (AMREX_SPACEDIM==1) + windir[2] = 1.0; + amrex::ignore_unused(dir); +#elif (AMREX_SPACEDIM==2) windir[2*dir] = 1.0; #else windir[dir] = 1.0; @@ -243,8 +255,10 @@ LaserParticleContainer::ContinuousInjection (const RealBox& injection_box) // Convert updated_position to Real* to use RealBox::contains(). #if (AMREX_SPACEDIM == 3) const Real* p_pos = m_updated_position.dataPtr(); -#else +#elif (AMREX_SPACEDIM == 2) const Real p_pos[2] = {m_updated_position[0], m_updated_position[2]}; +#else + const Real p_pos[1] = {m_updated_position[2]}; #endif if ( injection_box.contains(p_pos) ){ // Update laser_injection_box with current value @@ -278,6 +292,13 @@ LaserParticleContainer::UpdateContinuousInjectionPosition (Real dt) // which has 3 components, for both 2D and 3D simulations. m_updated_position[2*dir] -= WarpX::beta_boost * WarpX::boost_direction[2*dir] * PhysConst::c * dt; +#elif ( AMREX_SPACEDIM == 1 ) + // In 1D, dir=0 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for 1D, 2D, and 3D simulations. + m_updated_position[2] -= WarpX::beta_boost * + WarpX::boost_direction[2] * PhysConst::c * dt; + amrex::ignore_unused(dir); #endif } } @@ -317,12 +338,13 @@ LaserParticleContainer::InitData (int lev) m_position = m_updated_position; } +#if (AMREX_SPACEDIM >= 2) auto Transform = [&](int const i, int const j) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) return { m_position[0] + (S_X*(Real(i)+0.5_rt))*m_u_X[0] + (S_Y*(Real(j)+0.5_rt))*m_u_Y[0], m_position[1] + (S_X*(Real(i)+0.5_rt))*m_u_X[1] + (S_Y*(Real(j)+0.5_rt))*m_u_Y[1], m_position[2] + (S_X*(Real(i)+0.5_rt))*m_u_X[2] + (S_Y*(Real(j)+0.5_rt))*m_u_Y[2] }; -#else +#elif (AMREX_SPACEDIM == 2) amrex::ignore_unused(j); # if (defined WARPX_DIM_RZ) return { m_position[0] + (S_X*(Real(i)+0.5_rt)), @@ -335,18 +357,21 @@ LaserParticleContainer::InitData (int lev) # endif #endif }; +#endif // Given the "lab" frame coordinates, return the real coordinates in the laser plane coordinates auto InverseTransform = [&](const Vector<Real>& pos) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) return {m_u_X[0]*(pos[0]-m_position[0])+m_u_X[1]*(pos[1]-m_position[1])+m_u_X[2]*(pos[2]-m_position[2]), m_u_Y[0]*(pos[0]-m_position[0])+m_u_Y[1]*(pos[1]-m_position[1])+m_u_Y[2]*(pos[2]-m_position[2])}; -#else +#elif (AMREX_SPACEDIM == 2) # if (defined WARPX_DIM_RZ) return {pos[0]-m_position[0], 0.0_rt}; # else return {m_u_X[0]*(pos[0]-m_position[0])+m_u_X[2]*(pos[2]-m_position[2]), 0.0_rt}; # endif +#else + return {m_u_X[2]*(pos[2]-m_position[2]), 0.0_rt}; #endif }; @@ -374,11 +399,14 @@ LaserParticleContainer::InitData (int lev) compute_min_max(prob_hi[0], prob_lo[1], prob_hi[2]); compute_min_max(prob_lo[0], prob_hi[1], prob_hi[2]); compute_min_max(prob_hi[0], prob_hi[1], prob_hi[2]); -#else +#elif (AMREX_SPACEDIM == 2) compute_min_max(prob_lo[0], 0.0, prob_lo[1]); compute_min_max(prob_hi[0], 0.0, prob_lo[1]); compute_min_max(prob_lo[0], 0.0, prob_hi[1]); compute_min_max(prob_hi[0], 0.0, prob_hi[1]); +#else + compute_min_max(0.0, 0.0, prob_lo[0]); + compute_min_max(0.0, 0.0, prob_hi[0]); #endif } @@ -404,8 +432,10 @@ LaserParticleContainer::InitData (int lev) } } } -#else +#elif (AMREX_SPACEDIM == 2) BoxArray plane_ba { Box {IntVect(plane_lo[0],0), IntVect(plane_hi[0],0)} }; +#else + BoxArray plane_ba { Box {IntVect(0), IntVect(0)} }; #endif amrex::Vector<amrex::Real> particle_x, particle_y, particle_z, particle_w; @@ -419,11 +449,17 @@ LaserParticleContainer::InitData (int lev) const Box& bx = plane_ba[i]; for (IntVect cell = bx.smallEnd(); cell <= bx.bigEnd(); bx.next(cell)) { +#if (AMREX_SPACEDIM >= 2) const Vector<Real>& pos = Transform(cell[0], cell[1]); +#else + const Vector<Real>& pos = { 0.0_rt, 0.0_rt, m_position[2] }; +#endif #if (AMREX_SPACEDIM == 3) const Real* x = pos.dataPtr(); -#else +#elif (AMREX_SPACEDIM == 2) const Real x[2] = {pos[0], pos[2]}; +#else + const Real x[1] = {pos[2]}; #endif if (m_laser_injection_box.contains(x)) { @@ -583,6 +619,7 @@ LaserParticleContainer::Evolve (int lev, } } + if (rho && ! skip_deposition) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 1, 0, @@ -636,7 +673,7 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const Sy = std::min(std::min(dx[0]/(std::abs(m_u_Y[0])+eps), dx[1]/(std::abs(m_u_Y[1])+eps)), dx[2]/(std::abs(m_u_Y[2])+eps)); -#else +#elif (AMREX_SPACEDIM == 2) # if (defined WARPX_DIM_RZ) Sx = dx[0]; # else @@ -644,6 +681,10 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const dx[2]/(std::abs(m_u_X[2])+eps)); # endif Sy = 1.0; +#else + Sx = 1.0; + Sy = 1.0; + amrex::ignore_unused(eps); #endif } @@ -694,6 +735,7 @@ LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& p { const auto GetPosition = GetParticlePosition(pti); +#if (AMREX_SPACEDIM >= 2) Real tmp_u_X_0 = m_u_X[0]; Real tmp_u_X_2 = m_u_X[2]; Real tmp_position_0 = m_position[0]; @@ -705,6 +747,7 @@ LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& p Real tmp_u_Y_2 = m_u_Y[2]; Real tmp_position_1 = m_position[1]; #endif +#endif amrex::ParallelFor( np, @@ -725,6 +768,9 @@ LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& p tmp_u_X_0 * (x - tmp_position_0) + tmp_u_X_2 * (z - tmp_position_2); pplane_Yp[i] = 0.; +#elif (AMREX_SPACEDIM == 1) + pplane_Xp[i] = 0.; + pplane_Yp[i] = 0.; #endif } ); @@ -792,7 +838,9 @@ LaserParticleContainer::update_laser_particle (WarpXParIter& pti, // Push the the particle positions ParticleReal x, y, z; GetPosition(i, x, y, z); +#if !(defined WARPX_DIM_1D_Z) x += vx * dt; +#endif #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) y += vy * dt; #endif diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index b877708f9..70254aadd 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -357,6 +357,10 @@ MultiParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::use_fdtd_nci_corr==0, "ERROR: use_fdtd_nci_corr is not supported in RZ"); #endif +#ifdef WARPX_DIM_1D_Z + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::use_fdtd_nci_corr==0, + "ERROR: use_fdtd_nci_corr is not supported in 1D"); +#endif ParmParse pp_lasers("lasers"); pp_lasers.queryarr("names", lasers_names); @@ -1303,6 +1307,9 @@ MultiParticleContainer::doQEDSchwinger () #ifdef WARPX_DIM_RZ amrex::Abort("Schwinger process not implemented in rz geometry"); #endif +#ifdef WARPX_DIM_1D_Z + amrex::Abort("Schwinger process not implemented in 1D geometry"); +#endif // Get cell volume. In 2D the transverse size is // chosen by the user in the input file. diff --git a/Source/Particles/ParticleBoundaries_K.H b/Source/Particles/ParticleBoundaries_K.H index f33b29c17..48be4e65a 100644 --- a/Source/Particles/ParticleBoundaries_K.H +++ b/Source/Particles/ParticleBoundaries_K.H @@ -82,7 +82,10 @@ namespace ApplyParticleBoundaries { */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void - apply_boundaries (amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax, + apply_boundaries ( +#ifndef WARPX_DIM_1D_Z + amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax, +#endif #ifdef WARPX_DIM_3D amrex::ParticleReal& y, amrex::Real ymin, amrex::Real ymax, #endif @@ -96,10 +99,12 @@ namespace ApplyParticleBoundaries { bool change_sign_uy = false; bool change_sign_uz = false; +#ifndef WARPX_DIM_1D_Z apply_boundary(x, xmin, xmax, change_sign_ux, particle_lost, boundaries.xmin_bc, boundaries.xmax_bc, boundaries.reflection_model_xlo(-ux), boundaries.reflection_model_xhi(ux), engine); +#endif #ifdef WARPX_DIM_3D apply_boundary(y, ymin, ymax, change_sign_uy, particle_lost, boundaries.ymin_bc, boundaries.ymax_bc, diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp index 64efcefeb..d8856c297 100644 --- a/Source/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Particles/ParticleBoundaryBuffer.cpp @@ -69,12 +69,17 @@ ParticleBoundaryBuffer::ParticleBoundaryBuffer () for (int ispecies = 0; ispecies < numSpecies(); ++ispecies) { amrex::ParmParse pp_species(getSpeciesNames()[ispecies]); +#if AMREX_SPACEDIM == 1 + pp_species.query("save_particles_at_zlo", m_do_boundary_buffer[0][ispecies]); + pp_species.query("save_particles_at_zhi", m_do_boundary_buffer[1][ispecies]); +#elif AMREX_SPACEDIM == 2 pp_species.query("save_particles_at_xlo", m_do_boundary_buffer[0][ispecies]); pp_species.query("save_particles_at_xhi", m_do_boundary_buffer[1][ispecies]); -#if AMREX_SPACEDIM == 2 pp_species.query("save_particles_at_zlo", m_do_boundary_buffer[2][ispecies]); pp_species.query("save_particles_at_zhi", m_do_boundary_buffer[3][ispecies]); #else + pp_species.query("save_particles_at_xlo", m_do_boundary_buffer[0][ispecies]); + pp_species.query("save_particles_at_xhi", m_do_boundary_buffer[1][ispecies]); pp_species.query("save_particles_at_ylo", m_do_boundary_buffer[2][ispecies]); pp_species.query("save_particles_at_yhi", m_do_boundary_buffer[3][ispecies]); pp_species.query("save_particles_at_zlo", m_do_boundary_buffer[4][ispecies]); diff --git a/Source/Particles/ParticleCreation/SmartCreate.H b/Source/Particles/ParticleCreation/SmartCreate.H index 926168985..957efb2c4 100644 --- a/Source/Particles/ParticleCreation/SmartCreate.H +++ b/Source/Particles/ParticleCreation/SmartCreate.H @@ -46,13 +46,17 @@ struct SmartCreate const int cpu = 0, const int id = 0) const noexcept { - prt.m_aos[i_prt].pos(0) = x; #if (AMREX_SPACEDIM == 3) + prt.m_aos[i_prt].pos(0) = x; prt.m_aos[i_prt].pos(1) = y; prt.m_aos[i_prt].pos(2) = z; #elif (AMREX_SPACEDIM == 2) + prt.m_aos[i_prt].pos(0) = x; prt.m_aos[i_prt].pos(1) = z; amrex::ignore_unused(y); +#else + prt.m_aos[i_prt].pos(0) = z; + amrex::ignore_unused(x,y); #endif prt.m_aos[i_prt].cpu() = cpu; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c331b7a61..e34c4c28b 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -145,7 +145,7 @@ namespace pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; pos.y = lo_corner[1] + (iv[1]+r.y)*dx[1]; pos.z = lo_corner[2] + (iv[2]+r.z)*dx[2]; -#else +#elif (AMREX_SPACEDIM == 2) pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; pos.y = 0.0_rt; #if defined WARPX_DIM_XZ @@ -154,6 +154,10 @@ namespace // Note that for RZ, r.y will be theta pos.z = lo_corner[1] + (iv[1]+r.z)*dx[1]; #endif +#else + pos.x = 0.0_rt; + pos.y = 0.0_rt; + pos.z = lo_corner[0] + (iv[0]+r.z)*dx[0]; #endif return pos; } @@ -185,7 +189,9 @@ namespace ) noexcept { p.pos(0) = 0._rt; +#if (AMREX_SPACEDIM >= 2) p.pos(1) = 0._rt; +#endif #if (AMREX_SPACEDIM == 3) p.pos(2) = 0._rt; #endif @@ -283,7 +289,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If old particle positions should be saved add the needed components pp_species_name.query("save_previous_position", m_save_previous_position); if (m_save_previous_position) { +#if (AMREX_SPACEDIM >= 2) AddRealComp("prev_x"); +#endif #if (AMREX_SPACEDIM == 3) AddRealComp("prev_y"); #endif @@ -435,6 +443,11 @@ PhysicalParticleContainer::AddGaussianBeam ( const Real x = amrex::RandomNormal(x_m, x_rms); constexpr Real y = 0._prt; const Real z = amrex::RandomNormal(z_m, z_rms); +#elif (defined WARPX_DIM_1D_Z) + const Real weight = q_tot/(npart*charge*x_rms*y_rms); + constexpr Real x = 0._prt; + constexpr Real y = 0._prt; + const Real z = amrex::RandomNormal(z_m, z_rms); #endif if (plasma_injector->insideBounds(x, y, z) && std::abs( x - x_m ) < x_cut * x_rms && @@ -506,18 +519,20 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot, openPMD::ParticleSpecies ps = it.particles.begin()->second; auto const npart = ps["position"]["x"].getExtent()[0]; +#if !(defined WARPX_DIM_1D_Z) std::shared_ptr<ParticleReal> ptr_x = ps["position"]["x"].loadChunk<ParticleReal>(); double const position_unit_x = ps["position"]["x"].unitSI(); +#endif std::shared_ptr<ParticleReal> ptr_z = ps["position"]["z"].loadChunk<ParticleReal>(); double const position_unit_z = ps["position"]["z"].unitSI(); std::shared_ptr<ParticleReal> ptr_ux = ps["momentum"]["x"].loadChunk<ParticleReal>(); double const momentum_unit_x = ps["momentum"]["x"].unitSI(); std::shared_ptr<ParticleReal> ptr_uz = ps["momentum"]["z"].loadChunk<ParticleReal>(); double const momentum_unit_z = ps["momentum"]["z"].unitSI(); -# ifndef WARPX_DIM_XZ +# if !(defined(WARPX_DIM_XZ) || defined(WARPX_DIM_1D_Z)) std::shared_ptr<ParticleReal> ptr_y = ps["position"]["y"].loadChunk<ParticleReal>(); double const position_unit_y = ps["position"]["y"].unitSI(); -# endif +#endif std::shared_ptr<ParticleReal> ptr_uy = nullptr; double momentum_unit_y = 1.0; if (ps["momentum"].contains("y")) { @@ -547,13 +562,17 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot, } for (auto i = decltype(npart){0}; i<npart; ++i){ +#if !(defined WARPX_DIM_1D_Z) ParticleReal const x = ptr_x.get()[i]*position_unit_x; +#else + ParticleReal const x = 0.0_prt; +#endif ParticleReal const z = ptr_z.get()[i]*position_unit_z+z_shift; -# if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) ParticleReal const y = ptr_y.get()[i]*position_unit_y; -# else +#else ParticleReal const y = 0.0_prt; -# endif +#endif if (plasma_injector->insideBounds(x, y, z)) { ParticleReal const ux = ptr_ux.get()[i]*momentum_unit_x/PhysConst::m_e; ParticleReal const uz = ptr_uz.get()[i]*momentum_unit_z/PhysConst::m_e; @@ -711,6 +730,8 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) scale_fac = dx[0]*dx[1]*dx[2]/num_ppc; #elif AMREX_SPACEDIM==2 scale_fac = dx[0]*dx[1]/num_ppc; +#elif AMREX_SPACEDIM==1 + scale_fac = dx[0]/num_ppc; #endif } @@ -839,9 +860,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pcounts[index] = num_ppc; } } -#if (AMREX_SPACEDIM != 3) +#if (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); #endif +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(j,k); +#endif }); // Max number of new particles. All of them are created, @@ -954,7 +978,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) ); continue; } -#else +#elif (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi @@ -965,6 +989,17 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) ); continue; } +#else + amrex::ignore_unused(j,k); + if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi +#ifdef WARPX_QED + ,loc_has_quantum_sync, p_optical_depth_QSR + ,loc_has_breit_wheeler, p_optical_depth_BW +#endif + ); + continue; + } #endif // Save the x and y values to use in the insideBounds checks. @@ -1110,6 +1145,8 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #endif p.pos(0) = xb; p.pos(1) = pos.z; +#else + p.pos(0) = pos.z; #endif } }); @@ -1287,8 +1324,10 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) pcounts[index] = num_ppc_int; } } -#if (AMREX_SPACEDIM != 3) +#if (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); +#elif (AMREX_SPACEDIM == 1) + amrex::ignore_unused(j,k); #endif }); @@ -1402,12 +1441,18 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) p.id() = -1; continue; } -#else +#elif (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { p.id() = -1; continue; } +#else + amrex::ignore_unused(j,k); + if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { + p.id() = -1; + continue; + } #endif // Save the x and y values to use in the insideBounds checks. @@ -1491,6 +1536,8 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) #endif p.pos(0) = xb; p.pos(1) = pos.z; +#else + p.pos(0) = pos.z; #endif } }); @@ -1768,7 +1815,9 @@ PhysicalParticleContainer::applyNCIFilter ( const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::noz)}); +#elif (AMREX_SPACEDIM == 2) const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox), static_cast<int>(WarpX::noz)}); #else @@ -1836,7 +1885,13 @@ PhysicalParticleContainer::SplitParticles (int lev) long np_split; if(split_type==0) { - np_split = (AMREX_SPACEDIM == 3) ? 8 : 4; + if (AMREX_SPACEDIM == 3){ + np_split = 8; + } else if (AMREX_SPACEDIM == 2){ + np_split = 4; + } else { + np_split = 2; + } } else { np_split = 2*AMREX_SPACEDIM; } @@ -1876,7 +1931,20 @@ PhysicalParticleContainer::SplitParticles (int lev) // If particle is tagged, split it and put the // split particles in local arrays psplit_x etc. np_split_to_add += np_split; -#if (AMREX_SPACEDIM==2) +#if (AMREX_SPACEDIM==1) + // Split particle in two along z axis + // 2 particles in 1d, split_type doesn't matter? Discuss with Remi + for (int ishift = -1; ishift < 2; ishift +=2 ){ + // Add one particle with offset in z + psplit_x.push_back( xp ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + ishift*split_offset[2] ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + } +#elif (AMREX_SPACEDIM==2) if (split_type==0){ // Split particle in two along each diagonals // 4 particles in 2d @@ -2133,7 +2201,9 @@ PhysicalParticleContainer::GetParticleSlice ( WARPX_PROFILE("PhysicalParticleContainer::GetParticleSlice()"); // Assume that the boost in the positive z direction. -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + AMREX_ALWAYS_ASSERT(direction == 0); +#elif (AMREX_SPACEDIM == 2) AMREX_ALWAYS_ASSERT(direction == 1); #else AMREX_ALWAYS_ASSERT(direction == 2); @@ -2428,7 +2498,11 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, ParticleReal* y_old = nullptr; ParticleReal* z_old = nullptr; if (save_previous_position) { +#if (AMREX_SPACEDIM >= 2) x_old = pti.GetAttribs(particle_comps["prev_x"]).dataPtr(); +#else + amrex::ignore_unused(x_old); +#endif #if (AMREX_SPACEDIM == 3) y_old = pti.GetAttribs(particle_comps["prev_y"]).dataPtr(); #else @@ -2465,7 +2539,9 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, getPosition(ip, xp, yp, zp); if (save_previous_position) { +#if (AMREX_SPACEDIM >= 2) x_old[ip] = xp; +#endif #if (AMREX_SPACEDIM == 3) y_old[ip] = yp; #endif diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index b0f1257f5..4df1ef3c1 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -35,10 +35,14 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p, x = p.pos(0); y = p.pos(1); z = p.pos(2); -#else +#elif WARPX_DIM_XZ x = p.pos(0); y = amrex::ParticleReal(0.0); z = p.pos(1); +#else + x = amrex::ParticleReal(0.0); + y = amrex::ParticleReal(0.0); + z = p.pos(0); #endif } @@ -55,6 +59,9 @@ struct GetParticlePosition const RType* m_theta = nullptr; #elif (AMREX_SPACEDIM == 2) static constexpr RType m_y_default = RType(0.0); +#elif (AMREX_SPACEDIM == 1) + static constexpr RType m_x_default = RType(0.0); + static constexpr RType m_y_default = RType(0.0); #endif GetParticlePosition () = default; @@ -93,10 +100,14 @@ struct GetParticlePosition x = p.pos(0); y = p.pos(1); z = p.pos(2); -#else +#elif WARPX_DIM_XZ x = p.pos(0); y = m_y_default; z = p.pos(1); +#else + x = m_x_default; + y = m_y_default; + z = p.pos(0); #endif } @@ -117,10 +128,14 @@ struct GetParticlePosition x = p.pos(0); y = p.pos(1); z = p.pos(2); -#else +#elif WARPX_DIM_XZ x = p.pos(0); y = m_y_default; z = p.pos(1); +#else + x = m_x_default; + y = m_y_default; + z = p.pos(0); #endif } }; @@ -158,6 +173,9 @@ struct SetParticlePosition #if defined(WARPX_DIM_XZ) amrex::ignore_unused(y); #endif +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(x,y); +#endif #ifdef WARPX_DIM_RZ m_theta[i] = std::atan2(y, x); m_structs[i].pos(0) = std::sqrt(x*x + y*y); @@ -166,9 +184,11 @@ struct SetParticlePosition m_structs[i].pos(0) = x; m_structs[i].pos(1) = y; m_structs[i].pos(2) = z; -#else +#elif WARPX_DIM_XZ m_structs[i].pos(0) = x; m_structs[i].pos(1) = z; +#else + m_structs[i].pos(0) = z; #endif } @@ -182,6 +202,9 @@ struct SetParticlePosition #if defined(WARPX_DIM_XZ) amrex::ignore_unused(y); #endif +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(x,y); +#endif #ifdef WARPX_DIM_RZ m_structs[i].pos(0) = x; m_theta[i] = y; @@ -190,9 +213,11 @@ struct SetParticlePosition m_structs[i].pos(0) = x; m_structs[i].pos(1) = y; m_structs[i].pos(2) = z; -#else +#elif WARPX_DIM_XZ m_structs[i].pos(0) = x; m_structs[i].pos(1) = z; +#else + m_structs[i].pos(0) = z; #endif } }; diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 8ebb0b49a..39b12cd67 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -29,7 +29,11 @@ void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::Parti // Compute inverse Lorentz factor const amrex::Real inv_gamma = 1._rt/std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*inv_c2); // Update positions over one time step +#if (AMREX_SPACEDIM >= 2) x += ux * inv_gamma * dt; +#else + amrex::ignore_unused(x); +#endif #if (AMREX_SPACEDIM == 3) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * inv_gamma * dt; #else diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H index f52eb1c75..5e958c2c1 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -32,7 +32,11 @@ void UpdatePositionPhoton( const amrex::Real c_over_umod = (u_norm == 0._rt) ? 0._rt: PhysConst::c/u_norm; // Update positions over one time step +#if (AMREX_SPACEDIM >= 2) x += ux * c_over_umod * dt; +#else + amrex::ignore_unused(x); +#endif #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * c_over_umod * dt; #else diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 554d8fee4..f7cb27b28 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -109,12 +109,16 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) // The boundary conditions are read in in ReadBCParams but a child class // can allow these value to be overwritten if different boundary // conditions are desired for a specific species +#ifndef WARPX_DIM_1D_Z m_boundary_conditions.SetBoundsX(WarpX::particle_boundary_lo[0], WarpX::particle_boundary_hi[0]); +#endif #ifdef WARPX_DIM_3D m_boundary_conditions.SetBoundsY(WarpX::particle_boundary_lo[1], WarpX::particle_boundary_hi[1]); m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[2], WarpX::particle_boundary_hi[2]); -#else +#elif WARPX_DIM_XZ || WARPX_DIM_RZ m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[1], WarpX::particle_boundary_hi[1]); +#else + m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[0], WarpX::particle_boundary_hi[0]); #endif m_boundary_conditions.BuildReflectionModelParsers(); } @@ -212,6 +216,9 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, p.pos(0) = x[i]; #endif p.pos(1) = z[i]; +#else //AMREX_SPACEDIM == 1 + amrex::ignore_unused(x,y); + p.pos(0) = z[i]; #endif pinned_tile.push_back(p); @@ -325,7 +332,9 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, // deposit part of its current in a neighboring box. However, this should catch particles // traveling many cells away, for example with algorithms that allow for large time steps. -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::noz/2)); +#elif (AMREX_SPACEDIM == 2) const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::nox/2), static_cast<int>(WarpX::noz/2)); #elif (AMREX_SPACEDIM == 3) @@ -608,7 +617,9 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // are not trivial, this check might be too strict and we might need to relax it, as currently // done for the current deposition. -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::noz/2+1)); +#elif (AMREX_SPACEDIM == 2) const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::nox/2+1), static_cast<int>(WarpX::noz/2+1)); #elif (AMREX_SPACEDIM == 3) @@ -1132,8 +1143,10 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ { auto GetPosition = GetParticlePosition(pti); auto SetPosition = SetParticlePosition(pti); +#ifndef WARPX_DIM_1D_Z const Real xmin = Geom(lev).ProbLo(0); const Real xmax = Geom(lev).ProbHi(0); +#endif #ifdef WARPX_DIM_3D const Real ymin = Geom(lev).ProbLo(1); const Real ymax = Geom(lev).ProbHi(1); @@ -1163,7 +1176,10 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ // Note that for RZ, (x, y, z) is actually (r, theta, z). bool particle_lost = false; - ApplyParticleBoundaries::apply_boundaries(x, xmin, xmax, + ApplyParticleBoundaries::apply_boundaries( +#ifndef WARPX_DIM_1D_Z + x, xmin, xmax, +#endif #ifdef WARPX_DIM_3D y, ymin, ymax, #endif diff --git a/Source/Utils/CoarsenIO.cpp b/Source/Utils/CoarsenIO.cpp index 2cb97b0e1..e2ea419ad 100644 --- a/Source/Utils/CoarsenIO.cpp +++ b/Source/Utils/CoarsenIO.cpp @@ -41,24 +41,36 @@ CoarsenIO::Loop ( MultiFab& mf_dst, GpuArray<int,3> cr; // coarsening ratio sf[0] = stag_src[0]; +#if (AMREX_SPACEDIM == 1) + sf[1] = 0; +#else sf[1] = stag_src[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sf[2] = 0; #elif (AMREX_SPACEDIM == 3) sf[2] = stag_src[2]; #endif sc[0] = stag_dst[0]; +#if (AMREX_SPACEDIM == 1) + sc[1] = 0; +#else sc[1] = stag_dst[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sc[2] = 0; #elif (AMREX_SPACEDIM == 3) sc[2] = stag_dst[2]; #endif cr[0] = crse_ratio[0]; +#if (AMREX_SPACEDIM == 1) + cr[1] = 1; +#else cr[1] = crse_ratio[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) cr[2] = 1; #elif (AMREX_SPACEDIM == 3) cr[2] = crse_ratio[2]; diff --git a/Source/Utils/CoarsenMR.cpp b/Source/Utils/CoarsenMR.cpp index 80562d542..b701af088 100644 --- a/Source/Utils/CoarsenMR.cpp +++ b/Source/Utils/CoarsenMR.cpp @@ -30,24 +30,36 @@ CoarsenMR::Loop ( MultiFab& mf_dst, GpuArray<int,3> cr; // coarsening ratio sf[0] = stag_src[0]; +#if (AMREX_SPACEDIM == 1) + sf[1] = 0; +#else sf[1] = stag_src[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sf[2] = 0; #elif (AMREX_SPACEDIM == 3) sf[2] = stag_src[2]; #endif sc[0] = stag_dst[0]; +#if (AMREX_SPACEDIM == 1) + sc[1] = 0; +#else sc[1] = stag_dst[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sc[2] = 0; #elif (AMREX_SPACEDIM == 3) sc[2] = stag_dst[2]; #endif cr[0] = crse_ratio[0]; +#if (AMREX_SPACEDIM == 1) + cr[1] = 1; +#else cr[1] = crse_ratio[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) cr[2] = 1; #elif (AMREX_SPACEDIM == 3) cr[2] = crse_ratio[2]; diff --git a/Source/Utils/Interpolate_K.H b/Source/Utils/Interpolate_K.H index 9c8e0a2b4..f452f1acf 100644 --- a/Source/Utils/Interpolate_K.H +++ b/Source/Utils/Interpolate_K.H @@ -18,11 +18,16 @@ void interp (int j, int k, int l, Real const wx = static_cast<Real>(type[0]) * static_cast<amrex::Real>(j-jg*r_ratio) * rr; Real const owx = 1.0_rt-wx; +#if (AMREX_SPACEDIM >= 2) int const kg = amrex::coarsen(k,r_ratio); Real const wy = static_cast<Real>(type[1]) * static_cast<amrex::Real>(k-kg*r_ratio) * rr; Real const owy = 1.0_rt-wy; +#endif -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + fine(j,k,l) = owx * crse(jg ,0,0) + + wx * crse(jg+1,0,0); +#elif (AMREX_SPACEDIM == 2) fine(j,k,l) = owx * owy * crse(jg ,kg ,0) + owx * wy * crse(jg ,kg+1,0) + wx * owy * crse(jg+1,kg ,0) diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 600ad180e..3bb48d1a5 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -63,6 +63,12 @@ WarpX::UpdatePlasmaInjectionPosition (amrex::Real a_dt) // This needs to be converted in order to index `boost_direction` // which has 3 components, for both 2D and 3D simulations. WarpX::boost_direction[2*dir] * PhysConst::c * a_dt; +#elif ( AMREX_SPACEDIM == 1 ) + // In 1D, dir=0 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for 1D, 2D, and 3D simulations. + WarpX::boost_direction[2] * PhysConst::c * a_dt; + amrex::ignore_unused(dir); #endif } } @@ -365,13 +371,20 @@ WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // Compute x,y,z co-ordinates based on index type of mf +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0.0; + amrex::Real y = 0.0; + amrex::Real fac_z = (1.0 - mf_type[0]) * dx[0]*0.5; + amrex::Real z = i*dx[0] + real_box.lo(0) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1.0 - mf_type[0]) * dx[0]*0.5; amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0.0; amrex::Real fac_z = (1.0 - mf_type[1]) * dx[1]*0.5; amrex::Real z = j*dx[1] + real_box.lo(1) + fac_z; #else + amrex::Real fac_x = (1.0 - mf_type[0]) * dx[0]*0.5; + amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1.0 - mf_type[1]) * dx[1]*0.5; amrex::Real y = j*dx[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1.0 - mf_type[2]) * dx[2]*0.5; @@ -417,6 +430,11 @@ WarpX::ShiftGalileanBoundary () m_v_galilean[0]*time_shift, std::numeric_limits<amrex::Real>::quiet_NaN(), m_v_galilean[2]*time_shift }; +#elif (AMREX_SPACEDIM == 1) + m_galilean_shift = { + std::numeric_limits<Real>::quiet_NaN(), + std::numeric_limits<Real>::quiet_NaN(), + m_v_galilean[2]*time_shift }; #endif #if (AMREX_SPACEDIM == 3) @@ -430,8 +448,13 @@ WarpX::ShiftGalileanBoundary () new_hi[0] = current_hi[0] + m_galilean_shift[0]; new_lo[1] = current_lo[1] + m_galilean_shift[2]; new_hi[1] = current_hi[1] + m_galilean_shift[2]; - } - #endif + } +#elif (AMREX_SPACEDIM == 1) + { + new_lo[0] = current_lo[0] + m_galilean_shift[2]; + new_hi[0] = current_hi[0] + m_galilean_shift[2]; + } +#endif time_of_last_gal_shift = cur_time; ResetProbDomain(amrex::RealBox(new_lo, new_hi)); diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 06f0a7da5..d77a35b3a 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -170,8 +170,10 @@ void ConvertLabParamsToBoost() #if (AMREX_SPACEDIM == 3) Vector<int> dim_map {0, 1, 2}; -#else +#elif (AMREX_SPACEDIM == 2) Vector<int> dim_map {0, 2}; +#else + Vector<int> dim_map {2}; #endif for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) @@ -219,6 +221,8 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) amrex::Array<amrex::Real,3> galilean_shift = { 0., 0., 0., }; #elif (AMREX_SPACEDIM == 2) amrex::Array<amrex::Real,3> galilean_shift = { 0., std::numeric_limits<Real>::quiet_NaN(), 0., } ; +#elif (AMREX_SPACEDIM == 1) + amrex::Array<amrex::Real,3> galilean_shift = {std::numeric_limits<Real>::quiet_NaN(), std::numeric_limits<Real>::quiet_NaN(), 0., } ; #endif const amrex::Real zmin_box = WarpX::LowerCorner(bx, galilean_shift, lev)[2]; const amrex::Real zmax_box = WarpX::UpperCorner(bx, lev)[2]; @@ -226,8 +230,10 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) // Get box lower index in the z direction #if (AMREX_SPACEDIM==3) const int lo_ind = bx.loVect()[2]; -#else +#elif (AMREX_SPACEDIM==2) const int lo_ind = bx.loVect()[1]; +#else + const int lo_ind = bx.loVect()[0]; #endif // Check if box intersect with [zmin, zmax] if ( (zmax>zmin_box && zmin<=zmax_box) ){ @@ -237,8 +243,10 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept{ #if (AMREX_SPACEDIM==3) const Real z_gridpoint = zmin_box+(k-lo_ind)*dz; -#else +#elif (AMREX_SPACEDIM==2) const Real z_gridpoint = zmin_box+(j-lo_ind)*dz; +#else + const Real z_gridpoint = zmin_box+(i-lo_ind)*dz; #endif if ( (z_gridpoint >= zmin) && (z_gridpoint < zmax) ) { arr(i,j,k) = 0.; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index f41de8b58..5f9106388 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -659,7 +659,9 @@ WarpX::ReadParameters () Vector<int> parse_filter_npass_each_dir(AMREX_SPACEDIM,1); queryArrWithParser(pp_warpx, "filter_npass_each_dir", parse_filter_npass_each_dir, 0, AMREX_SPACEDIM); filter_npass_each_dir[0] = parse_filter_npass_each_dir[0]; +#if (AMREX_SPACEDIM >= 2) filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; +#endif #if (AMREX_SPACEDIM == 3) filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; #endif @@ -1432,7 +1434,9 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + amrex::RealVect dx({WarpX::CellSize(lev)[2]}); +#elif (AMREX_SPACEDIM == 2) amrex::RealVect dx = {WarpX::CellSize(lev)[0], WarpX::CellSize(lev)[2]}; #elif (AMREX_SPACEDIM == 3) amrex::RealVect dx = {WarpX::CellSize(lev)[0], WarpX::CellSize(lev)[1], WarpX::CellSize(lev)[2]}; @@ -1493,7 +1497,18 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm amrex::IntVect F_nodal_flag, G_nodal_flag; // Set nodal flags -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + // AMReX convention: x = missing dimension, y = missing dimension, z = only dimension + Ex_nodal_flag = IntVect(1); + Ey_nodal_flag = IntVect(1); + Ez_nodal_flag = IntVect(0); + Bx_nodal_flag = IntVect(0); + By_nodal_flag = IntVect(0); + Bz_nodal_flag = IntVect(1); + jx_nodal_flag = IntVect(1); + jy_nodal_flag = IntVect(1); + jz_nodal_flag = IntVect(0); +#elif (AMREX_SPACEDIM == 2) // AMReX convention: x = first dimension, y = missing dimension, z = second dimension Ex_nodal_flag = IntVect(0,1); Ey_nodal_flag = IntVect(1,1); @@ -2041,7 +2056,7 @@ WarpX::CellSize (int lev) #elif (AMREX_SPACEDIM == 2) return { dx[0], 1.0, dx[1] }; #else - static_assert(AMREX_SPACEDIM != 1, "1D is not supported"); + return { 1.0, 1.0, dx[0] }; #endif } @@ -2065,6 +2080,9 @@ WarpX::LowerCorner(const Box& bx, std::array<amrex::Real,3> galilean_shift, int #elif (AMREX_SPACEDIM == 2) return { xyzmin[0] + galilean_shift[0], std::numeric_limits<Real>::lowest(), xyzmin[1] + galilean_shift[2] }; + +#elif (AMREX_SPACEDIM == 1) + return { std::numeric_limits<Real>::lowest(), std::numeric_limits<Real>::lowest(), xyzmin[0] + galilean_shift[2] }; #endif } @@ -2077,6 +2095,8 @@ WarpX::UpperCorner(const Box& bx, int lev) return { xyzmax[0], xyzmax[1], xyzmax[2] }; #elif (AMREX_SPACEDIM == 2) return { xyzmax[0], std::numeric_limits<Real>::max(), xyzmax[1] }; +#elif (AMREX_SPACEDIM == 1) + return { std::numeric_limits<Real>::max(), std::numeric_limits<Real>::max(), xyzmax[0] }; #endif } @@ -2276,7 +2296,19 @@ WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_ma const auto hi = amrex::ubound( tbx ); Array4<int> msk = buffer_mask.array(); Array4<int const> gmsk = guard_mask.array(); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + int k = lo.z; + int j = lo.y; + for (int i = lo.x; i <= hi.x; ++i) { + setnull = false; + // If gmsk=0 for any neighbor within ng cells, current cell is in the buffer region + for (int ii = i-ng; ii <= i+ng; ++ii) { + if ( gmsk(ii,j,k) == 0 ) setnull = true; + } + if ( setnull ) msk(i,j,k) = 0; + else msk(i,j,k) = 1; + } +#elif (AMREX_SPACEDIM == 2) int k = lo.z; for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { diff --git a/cmake/WarpXFunctions.cmake b/cmake/WarpXFunctions.cmake index b5bbf0c9b..0d17d38ec 100644 --- a/cmake/WarpXFunctions.cmake +++ b/cmake/WarpXFunctions.cmake @@ -160,12 +160,10 @@ function(set_warpx_binary_name) endif() foreach(tgt IN LISTS _ALL_TARGETS) set_target_properties(${tgt} PROPERTIES OUTPUT_NAME "warpx") - if(WarpX_DIMS STREQUAL 3) - set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".3d") - elseif(WarpX_DIMS STREQUAL 2) - set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".2d") - elseif(WarpX_DIMS STREQUAL RZ) + if(WarpX_DIMS STREQUAL RZ) set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".RZ") + else() + set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".${WarpX_DIMS}d") endif() if(WarpX_MPI) @@ -226,12 +224,10 @@ function(set_warpx_binary_name) endif() if(WarpX_LIB) # alias to the latest build; this is the one expected by Python bindings - if(WarpX_DIMS STREQUAL 3) - set(lib_suffix "3d") - elseif(WarpX_DIMS STREQUAL 2) - set(lib_suffix "2d") - elseif(WarpX_DIMS STREQUAL RZ) + if(WarpX_DIMS STREQUAL RZ) set(lib_suffix "rz") + else() + set(lib_suffix "${WarpX_DIMS}d") endif() if(WIN32) set(mod_ext "dll") |