aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Example/Langmuir/inputs.multi.2d.rt86
-rw-r--r--Example/Langmuir/inputs.multi.rt12
-rwxr-xr-xExample/Langmuir/langmuir_multi_2d_analysis.py85
-rw-r--r--Regression/WarpX-tests.ini21
-rw-r--r--Source/WarpXFFT.cpp28
-rw-r--r--Source/WarpX_f.H13
-rw-r--r--Source/WarpX_fft.F905
-rw-r--r--tests/CurrentDeposition/GNUmakefile2
-rw-r--r--tests/FieldGather/GNUmakefile2
-rw-r--r--tests/FieldSolver/GNUmakefile2
-rw-r--r--tests/ParticlePusher/GNUmakefile2
11 files changed, 236 insertions, 22 deletions
diff --git a/Example/Langmuir/inputs.multi.2d.rt b/Example/Langmuir/inputs.multi.2d.rt
new file mode 100644
index 000000000..6f301a7c2
--- /dev/null
+++ b/Example/Langmuir/inputs.multi.2d.rt
@@ -0,0 +1,86 @@
+# Maximum number of time steps
+max_step = 80
+
+# number of grid points
+amr.n_cell = 128 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
+
+amr.plot_int = 40 # How often to write plotfiles. "<= 0" means no plotfiles.
+
+# Geometry
+geometry.coord_sys = 0 # 0: Cartesian
+geometry.is_periodic = 1 1 # Is periodic?
+geometry.prob_lo = -20.e-6 -20.e-6 # physical domain
+geometry.prob_hi = 20.e-6 20.e-6
+
+warpx.serialize_ics = 1
+
+# Verbosity
+warpx.verbose = 1
+
+# Algorithms
+algo.current_deposition = 3
+algo.charge_deposition = 0
+algo.field_gathering = 1
+algo.particle_pusher = 0
+
+# Interpolation
+interpolation.nox = 1
+interpolation.noy = 1
+interpolation.noz = 1
+
+# CFL
+warpx.cfl = 1.0
+
+# Parameters for the plasma wave
+constants.use_my_constants = 1
+constants.constant_names = epsilon kp k
+constants.constant_values = 0.01 376357.71524190728 314159.2653589793
+# 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.nspecies = 2
+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 2
+electrons.xmin = -20.e-6
+electrons.xmax = 20.e-6
+electrons.ymin = -20.e-6
+electrons.ymax = 20.e-6
+electrons.zmin = -20.e-6
+electrons.zmax = 20.e-6
+
+electrons.profile = constant
+electrons.density = 2.e24 # 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 2
+positrons.xmin = -20.e-6
+positrons.xmax = 20.e-6
+positrons.ymin = -20.e-6
+positrons.ymax = 20.e-6
+positrons.zmin = -20.e-6
+positrons.zmax = 20.e-6
+
+positrons.profile = constant
+positrons.density = 2.e24 # 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)"
diff --git a/Example/Langmuir/inputs.multi.rt b/Example/Langmuir/inputs.multi.rt
index edc5d3fa5..d2b648159 100644
--- a/Example/Langmuir/inputs.multi.rt
+++ b/Example/Langmuir/inputs.multi.rt
@@ -63,9 +63,9 @@ electrons.zmax = 20.e-6
electrons.profile = constant
electrons.density = 2.e24 # 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)"
+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
@@ -81,6 +81,6 @@ positrons.zmax = 20.e-6
positrons.profile = constant
positrons.density = 2.e24 # 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)"
+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)"
diff --git a/Example/Langmuir/langmuir_multi_2d_analysis.py b/Example/Langmuir/langmuir_multi_2d_analysis.py
new file mode 100755
index 000000000..1006cdd77
--- /dev/null
+++ b/Example/Langmuir/langmuir_multi_2d_analysis.py
@@ -0,0 +1,85 @@
+#! /usr/bin/env python
+
+# This is a script that analyses the simulation results from
+# the script `inputs.multi.rt`. This simulates a 3D periodic plasma wave.
+# The electric field in the simulation is given (in theory) by:
+# $$ E_x = \epsilon \,\frac{m_e c^2 k_x}{q_e}\sin(k_x x)\cos(k_y y)\cos(k_z z)\sin( \omega_p t)$$
+# $$ E_y = \epsilon \,\frac{m_e c^2 k_y}{q_e}\cos(k_x x)\sin(k_y y)\cos(k_z z)\sin( \omega_p t)$$
+# $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$
+import sys
+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
+
+# this will be the name of the plot file
+fn = sys.argv[1]
+
+# Parameters (these parameters must match the parameters in `inputs.multi.rt`)
+epsilon = 0.01
+n = 4.e24
+n_osc_x = 2
+n_osc_z = 2
+xmin = -20e-6; xmax = 20.e-6; Nx = 128
+zmin = -20e-6; zmax = 20.e-6; Nz = 128
+
+# Wave vector of the wave
+kx = 2.*np.pi*n_osc_x/(xmax-xmin)
+kz = 2.*np.pi*n_osc_z/(zmax-zmin)
+# Plasma frequency
+wp = np.sqrt((n*e**2)/(m_e*epsilon_0))
+
+k = {'Ex':kx, 'Ez':kz}
+cos = {'Ex': (0,1,1), 'Ez':(1,1,0)}
+
+def get_contribution( is_cos, k ):
+ du = (xmax-xmin)/Nx
+ u = xmin + du*( 0.5 + np.arange(Nx) )
+ 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]
+ x_contribution = get_contribution( cos_flag[0], kx )
+ z_contribution = get_contribution( cos_flag[2], kz )
+
+ E = amplitude * x_contribution[:, np.newaxis ] \
+ * z_contribution[np.newaxis, :]
+
+ return( E )
+
+# Read the file
+ds = yt.load(fn)
+t0 = ds.current_time.to_ndarray().mean()
+data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge,
+ dims=ds.domain_dimensions)
+
+# Check the validity of the fields
+overall_max_error = 0
+for field in ['Ex', 'Ez']:
+ E_sim = data[field].to_ndarray()[:,:,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))
+ overall_max_error = max( overall_max_error, max_error )
+
+# Plot the last field from the loop (Ez at iteration 40)
+plt.subplot2grid( (1,2), (0,0) )
+plt.imshow( E_sim )
+plt.colorbar()
+plt.title('Ez, last iteration\n(simulation)')
+plt.subplot2grid( (1,2), (0,1) )
+plt.imshow( E_th )
+plt.colorbar()
+plt.title('Ez, last iteration\n(theory)')
+plt.tight_layout()
+plt.savefig('langmuir_multi_2d_analysis.png')
+
+# Automatically check the validity
+assert overall_max_error < 0.02
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index b57e5c51d..c5e57f63c 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -147,6 +147,7 @@ numthreads = 2
compileTest = 0
doVis = 0
compareParticles = 1
+runtime_params = warpx.do_dynamic_scheduling=0
particleTypes = particle0 particle1
analysisRoutine = Example/Langmuir/langmuir_multi_analysis.py
analysisOutputImage = langmuir_multi_analysis.png
@@ -169,6 +170,24 @@ particleTypes = particle0 particle1
analysisRoutine = Example/Langmuir/langmuir_multi_analysis.py
analysisOutputImage = langmuir_multi_analysis.png
+[Langmuir_multi_2d_psatd]
+buildDir = .
+inputFile = Example/Langmuir/inputs.multi.2d.rt
+dim = 2
+addToCompileString = USE_PSATD=TRUE
+restartTest = 0
+useMPI = 1
+numprocs = 4
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 1
+runtime_params = psatd.fftw_plan_measure=0
+particleTypes = particle0 particle1
+analysisRoutine = Example/Langmuir/langmuir_multi_2d_analysis.py
+analysisOutputImage = langmuir_multi_2d_analysis.png
+
[LaserInjection]
buildDir = .
inputFile = Example/laser_injection/inputs.rt
@@ -250,7 +269,7 @@ dim = 3
restartTest = 0
useMPI = 1
numprocs = 2
-useOMP = 1
+useOMP = 0
numthreads = 0
compileTest = 0
doVis = 0
diff --git a/Source/WarpXFFT.cpp b/Source/WarpXFFT.cpp
index f394d5f64..810a1b29c 100644
--- a/Source/WarpXFFT.cpp
+++ b/Source/WarpXFFT.cpp
@@ -266,7 +266,6 @@ void
WarpX::FFTDomainDecompsition (int lev, BoxArray& ba_fft, DistributionMapping& dm_fft,
BoxArray& ba_valid, Box& domain_fft, const Box& domain)
{
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(AMREX_SPACEDIM == 3, "PSATD only works in 3D");
IntVect nguards_fft(AMREX_D_DECL(nox_fft/2,noy_fft/2,noz_fft/2));
@@ -281,14 +280,15 @@ WarpX::FFTDomainDecompsition (int lev, BoxArray& ba_fft, DistributionMapping& dm
// Ask FFTW to chop the current FFT sub-group domain in the z-direction
// and give a chunk to each MPI rank in the current sub-group.
int nz_fft, z0_fft;
- warpx_fft_domain_decomp(&nz_fft, &z0_fft, BL_TO_FORTRAN_BOX(domain_fft));
+
+ warpx_fft_domain_decomp(&nz_fft, &z0_fft, WARPX_TO_FORTRAN_BOX(domain_fft));
// Each MPI rank adds a box with its chunk of the FFT grid
// (given by the above decomposition) to the list `bx_fft`,
// then list is shared among all MPI ranks via AllGather
Vector<Box> bx_fft;
if (nz_fft > 0) {
Box b = domain_fft;
- b.setRange(2, z0_fft+domain_fft.smallEnd(2), nz_fft);
+ b.setRange(AMREX_SPACEDIM-1, z0_fft+domain_fft.smallEnd(AMREX_SPACEDIM-1), nz_fft);
bx_fft.push_back(b);
}
amrex::AllGatherBoxes(bx_fft);
@@ -397,17 +397,17 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */)
BL_PROFILE_VAR_START(blp_push_eb);
for (MFIter mfi(*Efield_fp_fft[lev][0]); mfi.isValid(); ++mfi)
{
- warpx_fft_push_eb(BL_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][0])[mfi]),
- BL_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][1])[mfi]),
- BL_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][2])[mfi]),
- BL_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][0])[mfi]),
- BL_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][1])[mfi]),
- BL_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][2])[mfi]),
- BL_TO_FORTRAN_ANYD((*current_fp_fft[lev][0])[mfi]),
- BL_TO_FORTRAN_ANYD((*current_fp_fft[lev][1])[mfi]),
- BL_TO_FORTRAN_ANYD((*current_fp_fft[lev][2])[mfi]),
- BL_TO_FORTRAN_ANYD((*rho_prev_fp_fft[lev])[mfi]),
- BL_TO_FORTRAN_ANYD((*rho_next_fp_fft[lev])[mfi]));
+ warpx_fft_push_eb(WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][0])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][1])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][2])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][0])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][1])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][2])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][0])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][1])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][2])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*rho_prev_fp_fft[lev])[mfi]),
+ WARPX_TO_FORTRAN_ANYD((*rho_next_fp_fft[lev])[mfi]));
}
BL_PROFILE_VAR_STOP(blp_push_eb);
diff --git a/Source/WarpX_f.H b/Source/WarpX_f.H
index 029c07377..377d8fd11 100644
--- a/Source/WarpX_f.H
+++ b/Source/WarpX_f.H
@@ -1,6 +1,19 @@
#include <AMReX_BLFort.H>
+#ifdef __cplusplus
+
+#if AMREX_SPACEDIM==2
+#define WARPX_ARLIM_ANYD(x) std::array<int,3>{(x)[0], 0, (x)[1]}.data()
+#else
+#define WARPX_ARLIM_ANYD(x) x
+#endif
+
+#define WARPX_TO_FORTRAN_BOX(x) WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect())
+#define WARPX_TO_FORTRAN_ANYD(x) (x).dataPtr(), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect())
+
+#endif
+
#if (AMREX_SPACEDIM == 3)
#define WRPX_COMPUTE_DIVB warpx_compute_divb_3d
diff --git a/Source/WarpX_fft.F90 b/Source/WarpX_fft.F90
index 406f3f90f..bdbdf7f3c 100644
--- a/Source/WarpX_fft.F90
+++ b/Source/WarpX_fft.F90
@@ -107,7 +107,6 @@ contains
type(c_ptr), intent(inout) :: fft_data(ndata)
real(c_double), intent(in) :: dx_wrpx(3), dt_wrpx
- integer :: iret
integer(idp) :: nopenmp
integer :: nx_padded
integer, dimension(3) :: shp
@@ -297,7 +296,11 @@ contains
rhoold => rhoold_wrpx
! Call the corresponding PICSAR function
+#if (BL_SPACEDIM == 3)
CALL push_psatd_ebfield_3d()
+#elif (BL_SPACEDIM == 2)
+ CALL push_psatd_ebfield_2d()
+#endif
ex => null()
ey => null()
diff --git a/tests/CurrentDeposition/GNUmakefile b/tests/CurrentDeposition/GNUmakefile
index bd400d019..05ac27a4b 100644
--- a/tests/CurrentDeposition/GNUmakefile
+++ b/tests/CurrentDeposition/GNUmakefile
@@ -17,6 +17,8 @@ include ./Make.package
include $(AMREX_HOME)/Src/Base/Make.package
include $(PICSAR_HOME)/src/Make.package
+DEFINES += -DWARPX
+
default: $(executable)
@echo SUCCESS
diff --git a/tests/FieldGather/GNUmakefile b/tests/FieldGather/GNUmakefile
index bd400d019..05ac27a4b 100644
--- a/tests/FieldGather/GNUmakefile
+++ b/tests/FieldGather/GNUmakefile
@@ -17,6 +17,8 @@ include ./Make.package
include $(AMREX_HOME)/Src/Base/Make.package
include $(PICSAR_HOME)/src/Make.package
+DEFINES += -DWARPX
+
default: $(executable)
@echo SUCCESS
diff --git a/tests/FieldSolver/GNUmakefile b/tests/FieldSolver/GNUmakefile
index bd400d019..05ac27a4b 100644
--- a/tests/FieldSolver/GNUmakefile
+++ b/tests/FieldSolver/GNUmakefile
@@ -17,6 +17,8 @@ include ./Make.package
include $(AMREX_HOME)/Src/Base/Make.package
include $(PICSAR_HOME)/src/Make.package
+DEFINES += -DWARPX
+
default: $(executable)
@echo SUCCESS
diff --git a/tests/ParticlePusher/GNUmakefile b/tests/ParticlePusher/GNUmakefile
index bd400d019..05ac27a4b 100644
--- a/tests/ParticlePusher/GNUmakefile
+++ b/tests/ParticlePusher/GNUmakefile
@@ -17,6 +17,8 @@ include ./Make.package
include $(AMREX_HOME)/Src/Base/Make.package
include $(PICSAR_HOME)/src/Make.package
+DEFINES += -DWARPX
+
default: $(executable)
@echo SUCCESS
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-13 03:30:06 +0000