diff options
21 files changed, 1618 insertions, 98 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 08804bb9b..3498b8203 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -540,6 +540,38 @@ Laser initialization none of the parameters below are used when ``<laser_name>.parse_field_function=1``. Even though ``<laser_name>.wavelength`` and ``<laser_name>.e_max`` should be included in the laser function, they still have to be specified as they are used for numerical purposes. + - ``"from_txye_file"``: the electric field of the laser is read from an external binary file + whose format is explained below. It requires to provide the name of the binary file + setting the additional parameter ``<laser_name>.txye_file_name`` (string). It accepts an + optional parameter ``<laser_name>.time_chunk_size`` (int). This allows to read only + time_chunk_size timesteps from the binary file. New timesteps are read as soon as they are needed. + The default value is automatically set to the number of timesteps contained in the binary file + (i.e. only one read is performed at the beginning of the simulation). + The external binary file should provide E(x,y,t) on a rectangular (but non necessarily uniform) + grid. The code performs a bi-linear (in 2D) or tri-linear (in 3D) interpolation to set the field + values. x,y,t are meant to be in S.I. units, while the field value is meant to be multiplied by + ``<laser_name>.e_max`` (i.e. in most cases the maximum of abs(E(x,y,t)) should be 1, + so that the maximum field intensity can be set straightforwardly with ``<laser_name>.e_max``). + The binary file has to respect the following format: + + * flag to indicate if the grid is uniform or not (1 byte, 0 means non-uniform, !=0 means uniform) + + * np, number of timesteps (uint32_t, must be >=2) + + * nx, number of points along x (uint32_t, must be >=2) + + * ny, number of points along y (uint32_t, must be 1 for 2D simulations and >=2 for 3D simulations) + + * timesteps (double[2] if grid is uniform, double[np] otherwise) + + * x_coords (double[2] if grid is uniform, double[nx] otherwise) + + * y_coords (double[1] if 2D, double[2] if 3D & uniform grid, double[ny] if 3D & non uniform grid) + + * field_data (double[nt * nx * ny], with nt being the slowest coordinate). + + A file at this format can be generated from Python, see an example at ``Examples/Modules/laser_injection_from_file`` + * ``<laser_name>.profile_t_peak`` (`float`; in seconds) The time at which the laser reaches its peak intensity, at the position @@ -637,7 +669,56 @@ Laser initialization ``mirror_z_width < dz/cell_size``, the upper bound of the mirror is increased so that it contains at least ``mirror_z_npoints``. -* ``warpx.E_external_grid`` & ``warpx.B_external_grid`` (list of `int`) optional (default `0. 0. 0.`) +* ``warpx.B_ext_grid_init_style`` (string) optional (default is "default") + This parameter determines the type of initialization for the external + magnetic field. The "default" style initializes the + external magnetic field (Bx,By,Bz) to (0.0, 0.0, 0.0). + The string can be set to "constant" if a constant magnetic field is + required to be set at initialization. If set to "constant", then an + additional parameter, namely, ``warpx.B_external_grid`` must be specified. + If set to ``parse_B_ext_grid_function``, then a mathematical expression can + be used to initialize the external magnetic field on the grid. It + requires additional parameters in the input file, namely, + ``warpx.Bx_external_grid_function(x,y,z)``, + ``warpx.By_external_grid_function(x,y,z)``, + ``warpx.Bz_external_grid_function(x,y,z)`` to initialize the external + magnetic field for each of the three components on the grid. + Constants required in the expression can be set using ``my_constants``. + For example, if ``warpx.Bx_external_grid_function(x,y,z)=Bo*x + delta*(y + z)`` + then the constants `Bo` and `delta` required in the above equation + can be set using ``my_constants.Bo=`` and ``my_constants.delta=`` in the + input file. For a two-dimensional simulation, it is assumed that the first dimension is `x` and the second dimension in `z`, and the value of `y` is set to zero. + Note that the current implementation of the parser for external B-field + does not work with RZ and the code will abort with an error message. + +* ``warpx.E_ext_grid_init_style`` (string) optional (default is "default") + This parameter determines the type of initialization for the external + electric field. The "default" style initializes the + external electric field (Ex,Ey,Ez) to (0.0, 0.0, 0.0). + The string can be set to "constant" if a constant electric field is + required to be set at initialization. If set to "constant", then an + additional parameter, namely, ``warpx.E_external_grid`` must be specified + in the input file. + If set to ``parse_E_ext_grid_function``, then a mathematical expression can + be used to initialize the external magnetic field on the grid. It + required additional parameters in the input file, namely, + ``warpx.Ex_external_grid_function(x,y,z)``, + ``warpx.Ey_external_grid_function(x,y,z)``, + ``warpx.Ez_externail_grid_function(x,y,z)`` to initialize the external + electric field for each of the three components on the grid. + Constants required in the expression can be set using ``my_constants``. + For example, if ``warpx.Ex_external_grid_function(x,y,z)=Eo*x + delta*(y + z)`` + then the constants `Bo` and `delta` required in the above equation + can be set using ``my_constants.Eo=`` and ``my_constants.delta=`` in the + input file. For a two-dimensional simulation, it is assumed that the first + dimension is `x` and the second dimension in `z`, + and the value of `y` is set to zero. + Note that the current implementation of the parser for external E-field + does not work with RZ and the code will abort with an error message. + +* ``warpx.E_external_grid`` & ``warpx.B_external_grid`` (list of `int`) + required when ``warpx.B_ext_grid_init_style="parse_B_ext_grid_function"`` + and when ``warpx.E_ext_grid_init_style="parse_E_ext_grid_function"``, respectively. External uniform and constant electrostatic and magnetostatic field added to the grid at initialization. Use with caution as these fields are used for the field solver. In particular, do not use any other boundary condition diff --git a/Examples/Modules/laser_injection_from_file/analysis.py b/Examples/Modules/laser_injection_from_file/analysis.py new file mode 100755 index 000000000..d1ef8a001 --- /dev/null +++ b/Examples/Modules/laser_injection_from_file/analysis.py @@ -0,0 +1,216 @@ +#!/usr/bin/python3 + +# This file is part of the WarpX automated test suite. It is used to test the +# injection of a laser pulse from an external binary file. +# +# - Generate an input binary file with a gaussian laser pulse. +# - Run the WarpX simulation for time T, when the pulse is fully injected +# - Compute the theory for laser envelope at time T +# - Compare theory and simulation, for both envelope and central frequency + +import yt ; yt.funcs.mylog.setLevel(50) +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import hilbert +import glob +import os + +#Maximum acceptable error for this test +relative_error_threshold = 0.06 + +#Physical parameters +um = 1.e-6 +fs = 1.e-15 +c = 299792458 + +#Parameters of the gaussian beam +wavelength = 1.*um +w0 = 6.*um +tt = 10.*fs +x_c = 0.*um +t_c = 20.*fs +foc_dist = 10*um +E_max = 1e12 +rot_angle = -np.pi/4.0 + +#Parameters of the tx grid +x_l = -12.0*um +x_r = 12.0*um +x_points = 480 +t_l = 0.0*fs +t_r = 40.0*fs +t_points = 400 +tcoords = np.linspace(t_l, t_r, t_points) +xcoords = np.linspace(x_l, x_r, x_points) + +def gauss(T,X,Y,opt): + """Compute the electric field for a Gaussian laser pulse. + This is used to write the binary input file. + """ + + k0 = 2.0*np.pi/wavelength + inv_tau2 = 1./tt/tt + osc_phase = k0*c*(T-t_c) + + diff_factor = 1.0 + 1.0j* foc_dist * 2/(k0*w0*w0) + inv_w_2 = 1.0/(w0*w0*diff_factor) + + pre_fact = np.exp(1.0j * osc_phase) + + if opt == '3d': + pre_fact = pre_fact/diff_factor + else: + pre_fact = pre_fact/np.sqrt(diff_factor) + + exp_arg = - (X*X + Y*Y)*inv_w_2 - inv_tau2 * (T-t_c)*(T-t_c) + + return np.real(pre_fact * np.exp(exp_arg)) + +# Function for the envelope +def gauss_env(T,XX,ZZ): + '''Function to compute the theory for the envelope + ''' + + X = np.cos(rot_angle)*XX + np.sin(rot_angle)*ZZ + Z = -np.sin(rot_angle)*XX + np.cos(rot_angle)*ZZ + + inv_tau2 = 1./tt/tt + inv_w_2 = 1.0/(w0*w0) + exp_arg = - (X*X)*inv_w_2 - inv_tau2 / c/c * (Z-T*c)*(Z-T*c) + return E_max * np.real(np.exp(exp_arg)) + +def write_file(fname, x, y, t, E): + """ For a given filename fname, space coordinates x and y, time coordinate t + and field E, write a WarpX-compatible input binary file containing the + profile of the laser pulse + """ + + with open(fname, 'wb') as file: + flag_unif = 0 + file.write(flag_unif.to_bytes(1, byteorder='little')) + file.write((len(t)).to_bytes(4, byteorder='little', signed=False)) + file.write((len(x)).to_bytes(4, byteorder='little', signed=False)) + file.write((len(y)).to_bytes(4, byteorder='little', signed=False)) + file.write(t.tobytes()) + file.write(x.tobytes()) + file.write(y.tobytes()) + file.write(E.tobytes()) + + +def write_file_unf(fname, x, y, t, E): + """ For a given filename fname, space coordinates x and y, time coordinate t + and field E, write a WarpX-compatible input binary file containing the + profile of the laser pulse. This function should be used in the case + of a uniform spatio-temporal mesh + """ + + with open(fname, 'wb') as file: + flag_unif = 1 + file.write(flag_unif.to_bytes(1, byteorder='little')) + file.write((len(t)).to_bytes(4, byteorder='little', signed=False)) + file.write((len(x)).to_bytes(4, byteorder='little', signed=False)) + file.write((len(y)).to_bytes(4, byteorder='little', signed=False)) + file.write(t[0].tobytes()) + file.write(t[-1].tobytes()) + file.write(x[0].tobytes()) + file.write(x[-1].tobytes()) + if len(y) == 1 : + file.write(y[0].tobytes()) + else : + file.write(y[0].tobytes()) + file.write(y[-1].tobytes()) + file.write(E.tobytes()) + +def create_gaussian_2d(): + T, X, Y = np.meshgrid(tcoords, xcoords, np.array([0.0]), indexing='ij') + E_t = gauss(T,X,Y,'2d') + write_file("gauss_2d.txye", xcoords, np.array([0.0]), tcoords, E_t) + write_file_unf("gauss_2d_unf.txye", xcoords, np.array([0.0]), tcoords, E_t) + + +def do_analysis(fname, compname, steps): + ds = yt.load(fname) + + dt = ds.current_time.to_value()/steps + + # Define 2D meshes + x = np.linspace( + ds.domain_left_edge[0], + ds.domain_right_edge[0], + ds.domain_dimensions[0]).v + z = np.linspace( + ds.domain_left_edge[ds.dimensionality-1], + ds.domain_right_edge[ds.dimensionality-1], + ds.domain_dimensions[ds.dimensionality-1]).v + X, Z = np.meshgrid(x, z, sparse=False, indexing='ij') + + # Compute the theory for envelope + env_theory = gauss_env(+t_c-ds.current_time.to_value(), X,Z)+gauss_env(-t_c+ds.current_time.to_value(), X,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) + F_laser = all_data_level_0['boxlib', 'Ey'].v.squeeze() + env = abs(hilbert(F_laser)) + extent = [ds.domain_left_edge[ds.dimensionality-1], ds.domain_right_edge[ds.dimensionality-1], + ds.domain_left_edge[0], ds.domain_right_edge[0] ] + + # Plot results + plt.figure(figsize=(8,6)) + plt.subplot(221) + plt.title('PIC field') + plt.imshow(F_laser, extent=extent) + plt.colorbar() + plt.subplot(222) + plt.title('PIC envelope') + plt.imshow(env, extent=extent) + plt.colorbar() + plt.subplot(223) + plt.title('Theory envelope') + plt.imshow(env_theory, extent=extent) + plt.colorbar() + plt.subplot(224) + plt.title('Difference') + plt.imshow(env-env_theory, extent=extent) + plt.colorbar() + plt.tight_layout() + plt.savefig(compname, bbox_inches='tight') + + relative_error_env = np.sum(np.abs(env-env_theory)) / np.sum(np.abs(env)) + print("Relative error envelope: ", relative_error_env) + assert(relative_error_env < relative_error_threshold) + + fft_F_laser = np.fft.fft2(F_laser) + + freq_rows = np.fft.fftfreq(F_laser.shape[0],dt) + freq_cols = np.fft.fftfreq(F_laser.shape[1],dt) + + pos_max = np.unravel_index(np.abs(fft_F_laser).argmax(), fft_F_laser.shape) + + freq = np.sqrt((freq_rows[pos_max[0]])**2 + (freq_cols[pos_max[1]]**2)) + exp_freq = c/wavelength + + relative_error_freq = np.abs(freq-exp_freq)/exp_freq + print("Relative error frequency: ", relative_error_freq) + assert(relative_error_freq < relative_error_threshold) + + + +def launch_analysis(executable): + create_gaussian_2d() + os.system("./" + executable + " inputs.2d_test_txye") + do_analysis("diags/plotfiles/plt00250/", "comp_unf.pdf", 250) + os.system("sed 's/gauss_2d_unf.txye/gauss_2d.txye/g' inputs.2d_test_txye > inputs.2d_test_txye_non_unf") + os.system("./" + executable + " inputs.2d_test_txye_non_unf") + do_analysis("diags/plotfiles/plt00250/", "comp_non_unf.pdf", 250) + + +def main() : + executables = glob.glob("main2d*") + if len(executables) == 1 : + launch_analysis(executables[0]) + else : + assert(False) + + +if __name__ == "__main__": + main() diff --git a/Examples/Modules/laser_injection_from_file/inputs.2d_test_txye b/Examples/Modules/laser_injection_from_file/inputs.2d_test_txye new file mode 100644 index 000000000..63ef9813a --- /dev/null +++ b/Examples/Modules/laser_injection_from_file/inputs.2d_test_txye @@ -0,0 +1,46 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 400 +amr.n_cell = 672 672 +amr.plot_int = 50 +amr.max_grid_size = 512 +amr.blocking_factor = 32 +amr.max_level = 0 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 1 1 # Is periodic? +geometry.prob_lo = -25.e-6 -25.0e-6 # physical domain +geometry.prob_hi = 25.e-6 25.e-6 +warpx.verbose = 1 +warpx.serialize_ics = 1 + +################################# +############ NUMERICS ########### +################################# +interpolation.nox = 3 +interpolation.noy = 3 +interpolation.noz = 3 +warpx.cfl = 0.98 +warpx.do_dynamic_scheduling = 0 +warpx.load_balance_int = -1 +warpx.use_filter = 0 +algo.maxwell_fdtd_solver = ckc + +################################# +############ PLASMA ############# +################################# +particles.nspecies = 0 + +################################# +############# LASER ############# +################################# +lasers.nlasers = 1 +lasers.names = txye_laser +txye_laser.position = 0. 0. 0. # This point is on the laser plane +txye_laser.direction = 1. 0. 1. # The plane normal direction +txye_laser.polarization = 0. 1. 0. # The main polarization vector +txye_laser.e_max = 1.e12 # Maximum amplitude of the laser field (in V/m) +txye_laser.wavelength = 1.0e-6 # The wavelength of the laser (in meters) +txye_laser.profile = from_txye_file +txye_laser.txye_file_name = "gauss_2d_unf.txye" +txye_laser.time_chunk_size = 50 diff --git a/Examples/Tests/collision/analysis_collision.py b/Examples/Tests/collision/analysis_collision.py index 25fa26e4b..06e8a3db0 100755 --- a/Examples/Tests/collision/analysis_collision.py +++ b/Examples/Tests/collision/analysis_collision.py @@ -19,7 +19,7 @@ import sys import yt import re import math -import statistics +import numpy from glob import glob tolerance = 0.001 @@ -53,8 +53,8 @@ for fn in fn_list: buf = temp.match(fn).groups() j = int(buf[1]) # compute error - vxe = statistics.mean(px[ 0:ne])/me/c - vxi = statistics.mean(px[ne:np])/mi/c + vxe = numpy.mean(px[ 0:ne])/me/c + vxi = numpy.mean(px[ne:np])/mi/c vxd = vxe - vxi fit = a*math.exp(b*j) error = error + abs(fit-vxd) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 7456c71a2..ed15fefa3 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -1007,6 +1007,20 @@ numthreads = 2 compileTest = 0 doVis = 0 +[LaserInjectionFromTXYEFile] +buildDir = . +inputFile = Examples/Modules/laser_injection_from_file/analysis.py +aux1File = Examples/Modules/laser_injection_from_file/inputs.2d_test_txye +customRunCmd = ./analysis.py +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 0 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 + [collision] buildDir = . inputFile = Examples/Tests/collision/inputs diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index f7c7e498f..5f75ed45a 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -3,6 +3,7 @@ #include <WarpXConst.H> #include <WarpX_f.H> #include <WarpX.H> +#include <WarpXUtil.H> #include <AMReX.H> @@ -168,34 +169,6 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) } } -namespace { -WarpXParser makeParser (std::string const& parse_function) -{ - WarpXParser parser(parse_function); - parser.registerVariables({"x","y","z"}); - - ParmParse pp("my_constants"); - std::set<std::string> symbols = parser.symbols(); - symbols.erase("x"); - symbols.erase("y"); - symbols.erase("z"); // after removing variables, we are left with constants - for (auto it = symbols.begin(); it != symbols.end(); ) { - Real v; - if (pp.query(it->c_str(), v)) { - parser.setConstant(*it, v); - it = symbols.erase(it); - } else { - ++it; - } - } - for (auto const& s : symbols) { // make sure there no unknown symbols - amrex::Abort("PlasmaInjector::makeParser: Unknown symbol "+s); - } - - return parser; -} -} - // Depending on injection type at runtime, initialize inj_rho // so that inj_rho->getDensity calls // InjectorPosition[Constant or Custom or etc.].getDensity. @@ -217,11 +190,7 @@ void PlasmaInjector::parseDensity (ParmParse& pp) // Construct InjectorDensity with InjectorDensityPredefined. inj_rho.reset(new InjectorDensity((InjectorDensityPredefined*)nullptr,species_name)); } else if (rho_prof_s == "parse_density_function") { - std::vector<std::string> f; - pp.getarr("density_function(x,y,z)", f); - for (auto const& s : f) { - str_density_function += s; - } + Store_parserString(pp, "density_function(x,y,z)", str_density_function); // Construct InjectorDensity with InjectorDensityParser. inj_rho.reset(new InjectorDensity((InjectorDensityParser*)nullptr, makeParser(str_density_function))); @@ -339,21 +308,12 @@ void PlasmaInjector::parseMomentum (ParmParse& pp) inj_mom.reset(new InjectorMomentum ((InjectorMomentumRadialExpansion*)nullptr, u_over_r)); } else if (mom_dist_s == "parse_momentum_function") { - std::vector<std::string> f; - pp.getarr("momentum_function_ux(x,y,z)", f); - for (auto const& s : f) { - str_momentum_function_ux += s; - } - f.clear(); - pp.getarr("momentum_function_uy(x,y,z)", f); - for (auto const& s : f) { - str_momentum_function_uy += s; - } - f.clear(); - pp.getarr("momentum_function_uz(x,y,z)", f); - for (auto const& s : f) { - str_momentum_function_uz += s; - } + Store_parserString(pp, "momentum_function_ux(x,y,z)", + str_momentum_function_ux); + Store_parserString(pp, "momentum_function_uy(x,y,z)", + str_momentum_function_uy); + Store_parserString(pp, "momentum_function_uz(x,y,z)", + str_momentum_function_uz); // Construct InjectorMomentum with InjectorMomentumParser. inj_mom.reset(new InjectorMomentum((InjectorMomentumParser*)nullptr, makeParser(str_momentum_function_ux), diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 29c9a8955..be29a1cbc 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1,4 +1,3 @@ - #include <WarpX.H> #include <WarpX_f.H> #include <BilinearFilter.H> @@ -10,6 +9,9 @@ #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> #endif +#include <GpuParser.H> +#include <WarpXUtil.H> + using namespace amrex; @@ -231,24 +233,161 @@ WarpX::PostRestart () mypc->PostRestart(); } + void WarpX::InitLevelData (int lev, Real time) { + + ParmParse pp("warpx"); + + // default values of E_external_grid and B_external_grid + // are used to set the E and B field when "constant" or + // "parser" is not explicitly used in the input. + pp.query("B_ext_grid_init_style", B_ext_grid_s); + std::transform(B_ext_grid_s.begin(), + B_ext_grid_s.end(), + B_ext_grid_s.begin(), + ::tolower); + + pp.query("E_ext_grid_init_style", E_ext_grid_s); + std::transform(E_ext_grid_s.begin(), + E_ext_grid_s.end(), + E_ext_grid_s.begin(), + ::tolower); + + // if the input string is "constant", the values for the + // external grid must be provided in the input. + if (B_ext_grid_s == "constant") + pp.getarr("B_external_grid", B_external_grid); + + // if the input string is "constant", the values for the + // external grid must be provided in the input. + if (E_ext_grid_s == "constant") + pp.getarr("E_external_grid", E_external_grid); + for (int i = 0; i < 3; ++i) { current_fp[lev][i]->setVal(0.0); - Efield_fp[lev][i]->setVal(E_external_grid[i]); - Bfield_fp[lev][i]->setVal(B_external_grid[i]); + if (lev > 0) + current_cp[lev][i]->setVal(0.0); + + if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") { + Bfield_fp[lev][i]->setVal(B_external_grid[i]); + if (lev > 0) { + Bfield_aux[lev][i]->setVal(B_external_grid[i]); + Bfield_cp[lev][i]->setVal(B_external_grid[i]); + } + } + if (E_ext_grid_s == "constant" || E_ext_grid_s == "default") { + Efield_fp[lev][i]->setVal(E_external_grid[i]); + if (lev > 0) { + Efield_aux[lev][i]->setVal(E_external_grid[i]); + Efield_cp[lev][i]->setVal(E_external_grid[i]); + } + } } - if (lev > 0) { - for (int i = 0; i < 3; ++i) { - Efield_aux[lev][i]->setVal(E_external_grid[i]); - Bfield_aux[lev][i]->setVal(B_external_grid[i]); + // if the input string for the B-field is "parse_b_ext_grid_function", + // then the analytical expression or function must be + // provided in the input file. + if (B_ext_grid_s == "parse_b_ext_grid_function") { - current_cp[lev][i]->setVal(0.0); - Efield_cp[lev][i]->setVal(E_external_grid[i]); - Bfield_cp[lev][i]->setVal(B_external_grid[i]); - } +#ifdef WARPX_DIM_RZ + amrex::Abort("E and B parser for external fields does not work with RZ -- TO DO"); +#endif + Store_parserString(pp, "Bx_external_grid_function(x,y,z)", + str_Bx_ext_grid_function); + Store_parserString(pp, "By_external_grid_function(x,y,z)", + str_By_ext_grid_function); + Store_parserString(pp, "Bz_external_grid_function(x,y,z)", + str_Bz_ext_grid_function); + + Bxfield_parser.reset(new ParserWrapper( + makeParser(str_Bx_ext_grid_function))); + Byfield_parser.reset(new ParserWrapper( + makeParser(str_By_ext_grid_function))); + Bzfield_parser.reset(new ParserWrapper( + makeParser(str_Bz_ext_grid_function))); + + // Initialize Bfield_fp with external function + InitializeExternalFieldsOnGridUsingParser(Bfield_fp[lev][0].get(), + Bfield_fp[lev][1].get(), + Bfield_fp[lev][2].get(), + Bxfield_parser.get(), + Byfield_parser.get(), + Bzfield_parser.get(), + Bx_nodal_flag, By_nodal_flag, + Bz_nodal_flag, lev); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Bfield_aux[lev][0].get(), + Bfield_aux[lev][1].get(), + Bfield_aux[lev][2].get(), + Bxfield_parser.get(), + Byfield_parser.get(), + Bzfield_parser.get(), + Bx_nodal_flag, By_nodal_flag, + Bz_nodal_flag, lev); + + InitializeExternalFieldsOnGridUsingParser(Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get(), + Bxfield_parser.get(), + Byfield_parser.get(), + Bzfield_parser.get(), + Bx_nodal_flag, By_nodal_flag, + Bz_nodal_flag, lev); + } + } + + // if the input string for the E-field is "parse_e_ext_grid_function", + // then the analytical expression or function must be + // provided in the input file. + if (E_ext_grid_s == "parse_e_ext_grid_function") { + +#ifdef WARPX_DIM_RZ + amrex::Abort("E and B parser for external fields does not work with RZ -- TO DO"); +#endif + Store_parserString(pp, "Ex_external_grid_function(x,y,z)", + str_Ex_ext_grid_function); + Store_parserString(pp, "Ey_external_grid_function(x,y,z)", + str_Ey_ext_grid_function); + Store_parserString(pp, "Ez_external_grid_function(x,y,z)", + str_Ez_ext_grid_function); + + Exfield_parser.reset(new ParserWrapper( + makeParser(str_Ex_ext_grid_function))); + Eyfield_parser.reset(new ParserWrapper( + makeParser(str_Ey_ext_grid_function))); + Ezfield_parser.reset(new ParserWrapper( + makeParser(str_Ez_ext_grid_function))); + + // Initialize Efield_fp with external function + InitializeExternalFieldsOnGridUsingParser(Efield_fp[lev][0].get(), + Efield_fp[lev][1].get(), + Efield_fp[lev][2].get(), + Exfield_parser.get(), + Eyfield_parser.get(), + Ezfield_parser.get(), + Ex_nodal_flag, Ey_nodal_flag, + Ez_nodal_flag, lev); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser(Efield_aux[lev][0].get(), + Efield_aux[lev][1].get(), + Efield_aux[lev][2].get(), + Exfield_parser.get(), + Eyfield_parser.get(), + Ezfield_parser.get(), + Ex_nodal_flag, Ey_nodal_flag, + Ez_nodal_flag, lev); + + InitializeExternalFieldsOnGridUsingParser(Efield_cp[lev][0].get(), + Efield_cp[lev][1].get(), + Efield_cp[lev][2].get(), + Exfield_parser.get(), + Eyfield_parser.get(), + Ezfield_parser.get(), + Ex_nodal_flag, Ey_nodal_flag, + Ez_nodal_flag, lev); + } } if (F_fp[lev]) { @@ -306,3 +445,100 @@ WarpX::InitLevelDataFFT (int lev, Real time) } #endif + +void +WarpX::InitializeExternalFieldsOnGridUsingParser ( + MultiFab *mfx, MultiFab *mfy, MultiFab *mfz, + ParserWrapper *xfield_parser, ParserWrapper *yfield_parser, + ParserWrapper *zfield_parser, IntVect x_nodal_flag, + IntVect y_nodal_flag, IntVect z_nodal_flag, + const int lev) +{ + + const auto dx_lev = geom[lev].CellSizeArray(); + const RealBox& real_box = geom[lev].ProbDomain(); + for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& tbx = mfi.tilebox(x_nodal_flag); + const Box& tby = mfi.tilebox(y_nodal_flag); + const Box& tbz = mfi.tilebox(z_nodal_flag); + + auto const& mfxfab = mfx->array(mfi); + auto const& mfyfab = mfy->array(mfi); + auto const& mfzfab = mfz->array(mfi); + + auto const& mfx_IndexType = (*mfx).ixType(); + auto const& mfy_IndexType = (*mfy).ixType(); + auto const& mfz_IndexType = (*mfz).ixType(); + + // Initialize IntVect based on the index type of multiFab + // 0 if cell-centered, 1 if node-centered. + IntVect mfx_type(AMREX_D_DECL(0,0,0)); + IntVect mfy_type(AMREX_D_DECL(0,0,0)); + IntVect mfz_type(AMREX_D_DECL(0,0,0)); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + mfx_type[idim] = mfx_IndexType.nodeCentered(idim); + mfy_type[idim] = mfy_IndexType.nodeCentered(idim); + mfz_type[idim] = mfz_IndexType.nodeCentered(idim); + } + + amrex::ParallelFor (tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + // Shift required in the x-, y-, or z- position + // depending on the index type of the multifab + Real fac_x = (1.0 - mfx_type[0]) * dx_lev[0]*0.5; + Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#else + Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mfx_type[2]) * dx_lev[2]*0.5; + Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif + // Initialize the x-component of the field. + mfxfab(i,j,k) = xfield_parser->getField(x,y,z); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real fac_x = (1.0 - mfy_type[0]) * dx_lev[0]*0.5; + Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==3) + Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mfx_type[2]) * dx_lev[2]*0.5; + Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif + // Initialize the y-component of the field. + mfyfab(i,j,k) = yfield_parser->getField(x,y,z); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real fac_x = (1.0 - mfz_type[0]) * dx_lev[0]*0.5; + Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==3) + Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mfz_type[2]) * dx_lev[2]*0.5; + Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif + // Initialize the z-component of the field. + mfzfab(i,j,k) = zfield_parser->getField(x,y,z); + }, + /* To allocate shared memory for the GPU threads. */ + /* But, for now only 3 doubles (x,y,z) are allocated. */ + amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3 + ); + + } + +} diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 1d9689979..ed9f5eda0 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -386,6 +386,9 @@ LaserParticleContainer::Evolve (int lev, t_lab = 1./WarpX::gamma_boost*t + WarpX::beta_boost*Z0_lab/PhysConst::c; } + // Update laser profile + m_up_laser_profile->update(t); + BL_ASSERT(OnSameGrids(lev,jx)); MultiFab* cost = WarpX::getCosts(lev); diff --git a/Source/Laser/LaserProfiles.H b/Source/Laser/LaserProfiles.H index e0ec0dc28..f97bf915e 100644 --- a/Source/Laser/LaserProfiles.H +++ b/Source/Laser/LaserProfiles.H @@ -5,6 +5,7 @@ #include <WarpXParser.H> #include <AMReX_ParmParse.H> #include <AMReX_Vector.H> +#include <AMReX_Gpu.H> #include <WarpXParser.H> @@ -13,6 +14,7 @@ #include <memory> #include <functional> #include <limits> +#include <utility> namespace WarpXLaserProfiles { @@ -31,8 +33,8 @@ struct CommonLaserParameters /** Abstract interface for laser profile classes * - * Each new laser profile should inherit this interface and implement two - * methods: init and fill_amplitude (described below). + * Each new laser profile should inherit this interface and implement three + * methods: init, update and fill_amplitude (described below). * * The declaration of a LaserProfile class should be placed in this file, * while the implementation of the methods should be in a dedicated file in @@ -60,6 +62,17 @@ public: const amrex::ParmParse& ppc, CommonLaserParameters params) = 0; + /** Update Laser Profile + * + * Some laser profiles might need to perform an "update" operation per + * time step. + * + * @param[in] t Current physical time in the simulation (seconds) + */ + virtual void + update ( + amrex::Real t) = 0; + /** Fill Electric Field Amplitude for each particle of the antenna. * * Xp, Yp and amplitude must be arrays with the same length @@ -76,7 +89,7 @@ public: amrex::Real const * AMREX_RESTRICT const Xp, amrex::Real const * AMREX_RESTRICT const Yp, amrex::Real t, - amrex::Real * AMREX_RESTRICT const amplitude) = 0; + amrex::Real * AMREX_RESTRICT const amplitude) const = 0; virtual ~ILaserProfile(){}; }; @@ -94,13 +107,17 @@ public: const amrex::ParmParse& ppc, CommonLaserParameters params) override final; + //No update needed + void + update (amrex::Real /*t */) override final {} + void fill_amplitude ( const int np, amrex::Real const * AMREX_RESTRICT const Xp, amrex::Real const * AMREX_RESTRICT const Yp, amrex::Real t, - amrex::Real * AMREX_RESTRICT const amplitude) override final; + amrex::Real * AMREX_RESTRICT const amplitude) const override final; private: struct { @@ -132,13 +149,17 @@ public: const amrex::ParmParse& ppc, CommonLaserParameters params) override final; + //No update needed + void + update (amrex::Real /*t */) override final {} + void fill_amplitude ( const int np, amrex::Real const * AMREX_RESTRICT const Xp, amrex::Real const * AMREX_RESTRICT const Yp, amrex::Real t, - amrex::Real * AMREX_RESTRICT const amplitude) override final; + amrex::Real * AMREX_RESTRICT const amplitude) const override final; private: struct { @@ -163,13 +184,17 @@ public: const amrex::ParmParse& ppc, CommonLaserParameters params) override final; + //No update needed + void + update (amrex::Real /*t */) override final {} + void fill_amplitude ( const int np, amrex::Real const * AMREX_RESTRICT const Xp, amrex::Real const * AMREX_RESTRICT const Yp, amrex::Real t, - amrex::Real * AMREX_RESTRICT const amplitude) override final; + amrex::Real * AMREX_RESTRICT const amplitude) const override final; private: struct{ @@ -180,6 +205,160 @@ private: }; /** + * Laser profile read from an 'TXYE' file + * The binary file must contain: + * - 3 unsigned integers (4 bytes): nt (points along t), nx (points along x) and ny (points along y) + * - nt*nx*ny doubles (8 bytes) in row major order : field amplitude + */ +class FromTXYEFileLaserProfile : public ILaserProfile +{ + +public: + void + init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& ppc, + CommonLaserParameters params) override final; + + /** \brief Reads new field data chunk from file if needed + * + * @param[in] t simulation time (seconds) + */ + void + update (amrex::Real t) override final; + + /** \brief compute field amplitude at particles' position for a laser beam + * loaded from an E(x,y,t) file. + * + * Both Xp and Yp are given in laser plane coordinate. + * For each particle with position Xp and Yp, this routine computes the + * amplitude of the laser electric field, stored in array amplitude. + * + * \param np: number of laser particles + * \param Xp: pointer to first component of positions of laser particles + * \param Yp: pointer to second component of positions of laser particles + * \param t: Current physical time + * \param amplitude: pointer to array of field amplitude. + */ + void + fill_amplitude ( + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) const override final; + + /** \brief Function to fill the amplitude in case of a uniform grid. + * This function cannot be private due to restrictions related to + * the use of extended __device__ lambda + * + * \param idx_t_left index of the last time coordinate < t + * \param np: number of laser particles + * \param Xp: pointer to first component of positions of laser particles + * \param Yp: pointer to second component of positions of laser particles + * \param t: Current physical time + * \param amplitude: pointer to array of field amplitude. + */ + void internal_fill_amplitude_uniform( + const int idx_t_left, + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) const; + + /** \brief Function to fill the amplitude in case of a non-uniform grid. + * This function cannot be private due to restrictions related to + * the use of extended __device__ lambda + * + * \param idx_t_left index of the last time coordinate < t + * \param np: number of laser particles + * \param Xp: pointer to first component of positions of laser particles + * \param Yp: pointer to second component of positions of laser particles + * \param t: Current physical time + * \param amplitude: pointer to array of field amplitude. + */ + void internal_fill_amplitude_nonuniform( + const int idx_t_left, + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) const; + +private: + /** \brief parse a field file in the binary 'txye' format (whose details are given below). + * + * A 'txye' file should be a binary file with the following format: + * -flag to indicate if the grid is uniform or not (1 byte, 0 means non-uniform, !=0 means uniform) + * -np, number of timesteps (uint32_t, must be >=2) + * -nx, number of points along x (uint32_t, must be >=2) + * -ny, number of points along y (uint32_t, must be 1 for 2D simulations and >=2 for 3D simulations) + * -timesteps (double[2] if grid is uniform, double[np] otherwise) + * -x_coords (double[2] if grid is uniform, double[nx] otherwise) + * -y_coords (double[1] if 2D, double[2] if 3D & uniform grid, double[ny] if 3D & non uniform grid) + * -field_data (double[nt * nx * ny], with nt being the slowest coordinate). + * The spatiotemporal grid must be rectangular. + * + * \param txye_file_name: name of the file to parse + */ + void parse_txye_file(std::string txye_file_name); + + /** \brief Finds left and right time indices corresponding to time t + * + * + * \param t: simulation time + */ + std::pair<int,int> find_left_right_time_indices(amrex::Real t) const; + + /** \brief Load field data within the temporal range [t_begin, t_end) + * + * Must be called after having parsed a data file with parse_txye_file. + * + * \param t_begin: left limit of the timestep range to read + * \param t_end: right limit of the timestep range to read (t_end is not read) + */ + void read_data_t_chuck(int t_begin, int t_end); + + /** + * \brief m_params contains all the internal parameters + * used by this laser profile + */ + struct{ + /** Name of the file containing the data */ + std::string txye_file_name; + /** Flag to store if the grid is uniform */ + bool is_grid_uniform = false; + /** Dimensions of E_data. nt, nx must be >=2. + * If DIM=3, ny must be >=2 as well. + * If DIM=2, ny must be 1 */ + int nt, nx, ny; + /** Vector of temporal coordinates. For a non-uniform grid, it contains + * all values of time. For a uniform grid, it contains only the start and stop + * times and intermediate times are obtained with nt */ + amrex::Gpu::ManagedVector<amrex::Real> t_coords; + /** Vector or x coordinates. For a non-uniform grid, it contains all values + * of space dimension x. For a uniform grid, it contains only the min and max + * x coordinates, and intermediate positions are obtained with nx */ + amrex::Gpu::ManagedVector<amrex::Real> x_coords; + /** Vector or y coordinates. For a non-uniform grid, it contains all values + * of space dimension y. For a uniform grid, it contains only the min and max + * y coordinates, and intermediate positions are obtained with ny */ + amrex::Gpu::ManagedVector<amrex::Real> y_coords; + /** Size of the timestep range to load */ + int time_chunk_size; + /** Index of the first timestep in memory */ + int first_time_index; + /** Index of the last timestep in memory */ + int last_time_index; + /** Field data */ + amrex::Gpu::ManagedVector<amrex::Real> E_data; + } m_params; + + CommonLaserParameters m_common_params; +}; + +/** * Maps laser profile names to lambdas returing unique pointers * to the corresponding laser profile objects. */ @@ -195,7 +374,9 @@ laser_profiles_dictionary = {"harris", [] () {return std::make_unique<HarrisLaserProfile>();} }, {"parse_field_function", - [] () {return std::make_unique<FieldFunctionLaserProfile>();} } + [] () {return std::make_unique<FieldFunctionLaserProfile>();} }, + {"from_txye_file", + [] () {return std::make_unique<FromTXYEFileLaserProfile>();} } }; } //WarpXLaserProfiles diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp index 3c9d7379a..d34bc6aba 100644 --- a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp +++ b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp @@ -37,7 +37,7 @@ FieldFunctionLaserProfile::init ( void FieldFunctionLaserProfile::fill_amplitude ( const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, - Real t, Real * AMREX_RESTRICT const amplitude) + Real t, Real * AMREX_RESTRICT const amplitude) const { for (int i = 0; i < np; ++i) { amplitude[i] = m_parser.eval(Xp[i], Yp[i], t); diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp new file mode 100644 index 000000000..8f44c2d5f --- /dev/null +++ b/Source/Laser/LaserProfilesImpl/LaserProfileFromTXYEFile.cpp @@ -0,0 +1,479 @@ +#include <LaserProfiles.H> + +#include <WarpX_Complex.H> +#include <WarpXConst.H> +#include <WarpXUtil.H> + +#include <AMReX_Print.H> +#include <AMReX_ParallelDescriptor.H> + +#include <limits> +#include <iostream> +#include <fstream> +#include <cstdint> +#include <algorithm> + +using namespace amrex; +using namespace WarpXLaserProfiles; + +void +FromTXYEFileLaserProfile::init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& /* ppc */, + CommonLaserParameters params) +{ + if (!std::numeric_limits< double >::is_iec559) + { + Print() << R"(Warning: double does not comply with IEEE 754: bad + things will happen parsing the X, Y and T profiles for the laser!)"; + } + + // Parse the TXYE file + ppl.get("txye_file_name", m_params.txye_file_name); + if(m_params.txye_file_name.empty()) + { + Abort("txye_file_name must be provided for txye_file laser profile!"); + } + parse_txye_file(m_params.txye_file_name); + + //Set time_chunk_size + m_params.time_chunk_size = m_params.nt; + int temp = 1; + if(ppl.query("time_chunk_size", temp)){ + m_params.time_chunk_size = min( + temp, m_params.time_chunk_size); + } + if(m_params.time_chunk_size < 2){ + Abort("Error! time_chunk_size must be >= 2!"); + } + + //Allocate memory for E_data Vector + const int data_size = m_params.time_chunk_size* + m_params.nx*m_params.ny; + m_params.E_data = Gpu::ManagedVector<amrex::Real>(data_size); + + //Read first time chunck + read_data_t_chuck(0, m_params.time_chunk_size); + + //Copy common params + m_common_params = params; +} + +void +FromTXYEFileLaserProfile::update (amrex::Real t) +{ + if(t >= m_params.t_coords.back()) + return; + + const auto idx_times = find_left_right_time_indices(t); + const auto idx_t_left = idx_times.first; + const auto idx_t_right = idx_times.second; + + //Load data chunck if needed + if(idx_t_right > m_params.last_time_index){ + read_data_t_chuck(idx_t_left, idx_t_left+m_params.time_chunk_size); + } +} + +void +FromTXYEFileLaserProfile::fill_amplitude ( + const int np, + Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, + Real t, Real * AMREX_RESTRICT const amplitude) const +{ + //Amplitude is 0 if time is out of range + if(t < m_params.t_coords.front() || t > m_params.t_coords.back()){ + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE (int i) { + amplitude[i] = 0.0_rt;}); + return; + } + + //Find left and right time indices + int idx_t_left, idx_t_right; + std::tie(idx_t_left, idx_t_right) = find_left_right_time_indices(t); + + if(idx_t_left < m_params.first_time_index){ + Abort("Something bad has happened with the simulation time"); + } + + if(m_params.is_grid_uniform){ + internal_fill_amplitude_uniform( + idx_t_left, np, Xp, Yp, t, amplitude); + } + else{ + internal_fill_amplitude_nonuniform( + idx_t_left, np, Xp, Yp, t, amplitude); + } +} + +void +FromTXYEFileLaserProfile::parse_txye_file(std::string txye_file_name) +{ + if(ParallelDescriptor::IOProcessor()){ + std::ifstream inp(txye_file_name, std::ios::binary); + if(!inp) Abort("Failed to open txye file"); + + //Uniform grid flag + char flag; + inp.read(&flag, 1); + if(!inp) Abort("Failed to read sizes from txye file"); + m_params.is_grid_uniform=flag; + + //Grid points along t, x and y + auto const three_uint32_size = sizeof(uint32_t)*3; + char buf[three_uint32_size]; + inp.read(buf, three_uint32_size); + if(!inp) Abort("Failed to read sizes from txye file"); + m_params.nt = reinterpret_cast<uint32_t*>(buf)[0]; + m_params.nx = reinterpret_cast<uint32_t*>(buf)[1]; + m_params.ny = reinterpret_cast<uint32_t*>(buf)[2]; + if(m_params.nt <= 1) Abort("nt in txye file must be >=2"); + if(m_params.nx <= 1) Abort("nx in txye file must be >=2"); +#if (AMREX_SPACEDIM == 3) + if(m_params.ny <= 1) Abort("ny in txye file must be >=2 in 3D"); +#elif(AMREX_SPACEDIM == 2) + if(m_params.ny != 1) Abort("ny in txye file must be 1 in 2D"); +#endif + + //Coordinates + Vector<double> buf_t, buf_x, buf_y; + if(m_params.is_grid_uniform){ + buf_t.resize(2); + buf_x.resize(2); +#if (AMREX_SPACEDIM == 3) + buf_y.resize(2); +#elif(AMREX_SPACEDIM == 2) + buf_y.resize(1); +#endif + } + else{ + buf_t.resize(m_params.nt); + buf_x.resize(m_params.nx); + buf_y.resize(m_params.ny); + } + inp.read(reinterpret_cast<char*>(buf_t.dataPtr()), + buf_t.size()*sizeof(double)); + if(!inp) + Abort("Failed to read coords from txye file"); + if (!std::is_sorted(buf_t.begin(), buf_t.end())) + Abort("Coordinates are not sorted in txye file"); + inp.read(reinterpret_cast<char*>(buf_x.dataPtr()), + buf_x.size()*sizeof(double)); + if(!inp) + Abort("Failed to read coords from txye file"); + if (!std::is_sorted(buf_x.begin(), buf_x.end())) + Abort("Coordinates are not sorted in txye file"); + inp.read(reinterpret_cast<char*>(buf_y.dataPtr()), + buf_y.size()*sizeof(double)); + if(!inp) + Abort("Failed to read coords from txye file"); + if (!std::is_sorted(buf_y.begin(), buf_y.end())) + Abort("Coordinates are not sorted in txye file"); + m_params.t_coords = Gpu::ManagedVector<amrex::Real>(buf_t.size()); + m_params.x_coords = Gpu::ManagedVector<amrex::Real>(buf_x.size()); + m_params.y_coords = Gpu::ManagedVector<amrex::Real>(buf_y.size()); + // Convert from double to amrex::Real + std::transform(buf_t.begin(), buf_t.end(), m_params.t_coords.begin(), + [](auto x) {return static_cast<amrex::Real>(x);} ); + std::transform(buf_x.begin(), buf_x.end(), m_params.x_coords.begin(), + [](auto x) {return static_cast<amrex::Real>(x);} ); + std::transform(buf_y.begin(), buf_y.end(), m_params.y_coords.begin(), + [](auto x) {return static_cast<amrex::Real>(x);} ); + } + + //Broadcast grid uniformity + char is_grid_uniform = m_params.is_grid_uniform; + ParallelDescriptor::Bcast(&is_grid_uniform, 1, + ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Barrier(); + m_params.is_grid_uniform = is_grid_uniform; + + //Broadcast grid size and coordinate sizes + //When a non-uniform grid is used, nt, nx and ny are identical + //to t_coords.size(), x_coords.size() and y_coords.size(). + //When a uniform grid is used, nt,nx and ny store the number of points + //used for the mesh, while t_coords, x_coords and y_coords store the + //extrems in each direaction. Thus t_coords and x_coords in this case + //have size 2 and y_coords has size 1 in 2D and size 2 in 3D. + int t_sizes[6] = {m_params.nt, m_params.nx, m_params.ny, + static_cast<int>(m_params.t_coords.size()), + static_cast<int>(m_params.x_coords.size()), + static_cast<int>(m_params.y_coords.size())}; + ParallelDescriptor::Bcast(t_sizes, 6, + ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Barrier(); + m_params.nt = t_sizes[0]; m_params.nx = t_sizes[1]; m_params.ny = t_sizes[2]; + + //Broadcast coordinates + if(!ParallelDescriptor::IOProcessor()){ + m_params.t_coords = Gpu::ManagedVector<amrex::Real>(t_sizes[3]); + m_params.x_coords = Gpu::ManagedVector<amrex::Real>(t_sizes[4]); + m_params.y_coords = Gpu::ManagedVector<amrex::Real>(t_sizes[5]); + } + ParallelDescriptor::Bcast(m_params.t_coords.dataPtr(), + m_params.t_coords.size(), ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Bcast(m_params.x_coords.dataPtr(), + m_params.x_coords.size(), ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Bcast(m_params.y_coords.dataPtr(), + m_params.y_coords.size(), ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Barrier(); +} + +std::pair<int,int> +FromTXYEFileLaserProfile::find_left_right_time_indices(amrex::Real t) const +{ + int idx_t_right; + if(m_params.is_grid_uniform){ + const auto t_min = m_params.t_coords.front(); + const auto t_max = m_params.t_coords.back(); + const auto temp_idx_t_right = static_cast<int>( + ceil( (m_params.nt-1)*(t-t_min)/(t_max-t_min))); + idx_t_right = max(min(temp_idx_t_right, m_params.nt-1),1); + } + else{ + idx_t_right = std::distance(m_params.t_coords.begin(), + std::upper_bound(m_params.t_coords.begin(), + m_params.t_coords.end(), t)); + } + return std::make_pair(idx_t_right-1, idx_t_right); +} + +void +FromTXYEFileLaserProfile::read_data_t_chuck(int t_begin, int t_end) +{ + amrex::Print() << + "Reading [" << t_begin << ", " << t_end << + ") data chunk from " << m_params.txye_file_name << "\n"; + + //Indices of the first and last timestep to read + auto i_first = max(0, t_begin); + auto i_last = min(t_end-1, m_params.nt-1); + if(i_last-i_first+1 > m_params.E_data.size()) + Abort("Data chunk to read from file is too large"); + + if(ParallelDescriptor::IOProcessor()){ + //Read data chunk + std::ifstream inp(m_params.txye_file_name, std::ios::binary); + if(!inp) Abort("Failed to open txye file"); + auto skip_amount = 1 + + 3*sizeof(uint32_t) + + m_params.t_coords.size()*sizeof(double) + + m_params.x_coords.size()*sizeof(double) + + m_params.y_coords.size()*sizeof(double) + + sizeof(double)*t_begin*m_params.nx*m_params.ny; + inp.ignore(skip_amount); + if(!inp) Abort("Failed to read field data from txye file"); + const int read_size = (i_last - i_first + 1)* + m_params.nx*m_params.ny; + Vector<double> buf_e(read_size); + inp.read(reinterpret_cast<char*>(buf_e.dataPtr()), read_size*sizeof(double)); + if(!inp) Abort("Failed to read field data from txye file"); + std::transform(buf_e.begin(), buf_e.end(), m_params.E_data.begin(), + [](auto x) {return static_cast<amrex::Real>(x);} ); + } + + //Broadcast E_data + ParallelDescriptor::Bcast(m_params.E_data.dataPtr(), + m_params.E_data.size(), ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Barrier(); + + //Update first and last indices + m_params.first_time_index = i_first; + m_params.last_time_index = i_last; +} + +void +FromTXYEFileLaserProfile::internal_fill_amplitude_uniform( + const int idx_t_left, + const int np, + Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, + Real t, Real * AMREX_RESTRICT const amplitude) const +{ + // Copy member variables to tmp copies + // and get pointers to underlying data for GPU. + const auto tmp_e_max = m_common_params.e_max; + const auto tmp_x_min = m_params.x_coords.front(); + const auto tmp_x_max = m_params.x_coords.back(); + const auto tmp_y_min = m_params.y_coords.front(); + const auto tmp_y_max = m_params.y_coords.back(); + const auto tmp_nx = m_params.nx; +#if (AMREX_SPACEDIM == 3) + const auto tmp_ny = m_params.ny; +#endif + const auto p_E_data = m_params.E_data.dataPtr(); + const auto tmp_idx_first_time = m_params.first_time_index; + const int idx_t_right = idx_t_left+1; + const auto t_left = idx_t_left* + (m_params.t_coords.back()-m_params.t_coords.front())/(m_params.nt-1) + + m_params.t_coords.front(); + const auto t_right = idx_t_right* + (m_params.t_coords.back()-m_params.t_coords.front())/(m_params.nt-1) + + m_params.t_coords.front(); + + // Loop through the macroparticle to calculate the proper amplitude + amrex::ParallelFor( + np, + [=] AMREX_GPU_DEVICE (int i) { + //Amplitude is zero if we are out of bounds + if (Xp[i] <= tmp_x_min || Xp[i] >= tmp_x_max){ + amplitude[i] = 0.0_rt; + return; + } +#if (AMREX_SPACEDIM == 3) + if (Yp[i] <= tmp_y_min || Yp[i] >= tmp_y_max){ + amplitude[i] = 0.0_rt; + return; + } +#endif + //Find indices and coordinates along x + const int temp_idx_x_right = static_cast<int>( + ceil((tmp_nx-1)*(Xp[i]- tmp_x_min)/(tmp_x_max-tmp_x_min))); + const int idx_x_right = + max(min(temp_idx_x_right,tmp_nx-1),static_cast<int>(1)); + const int idx_x_left = idx_x_right - 1; + const auto x_0 = + idx_x_left*(tmp_x_max-tmp_x_min)/(tmp_nx-1) + tmp_x_min; + const auto x_1 = + idx_x_right*(tmp_x_max-tmp_x_min)/(tmp_nx-1) + tmp_x_min; + +#if (AMREX_SPACEDIM == 2) + //Interpolate amplitude + const auto idx = [=](int i, int j){ + return (i-tmp_idx_first_time) * tmp_nx + j; + }; + amplitude[i] = WarpXUtilAlgo::bilinear_interp( + t_left, t_right, + x_0, x_1, + p_E_data[idx(idx_t_left, idx_x_left)], + p_E_data[idx(idx_t_left, idx_x_right)], + p_E_data[idx(idx_t_right, idx_x_left)], + p_E_data[idx(idx_t_right, idx_x_right)], + t, Xp[i])*tmp_e_max; + +#elif (AMREX_SPACEDIM == 3) + //Find indices and coordinates along y + const int temp_idx_y_right = static_cast<int>( + ceil((tmp_ny-1)*(Yp[i]- tmp_y_min)/(tmp_y_max-tmp_y_min))); + const int idx_y_right = + max(min(temp_idx_y_right,tmp_ny-1),static_cast<int>(1)); + const int idx_y_left = idx_y_right - 1; + const auto y_0 = + idx_y_left*(tmp_y_max-tmp_y_min)/(tmp_ny-1) + tmp_y_min; + const auto y_1 = + idx_y_right*(tmp_y_max-tmp_y_min)/(tmp_ny-1) + tmp_y_min; + + //Interpolate amplitude + const auto idx = [=](int i, int j, int k){ + return + (i-tmp_idx_first_time)*tmp_nx*tmp_ny+ + j*tmp_ny + k; + }; + amplitude[i] = WarpXUtilAlgo::trilinear_interp( + t_left, t_right, + x_0, x_1, + y_0, y_1, + p_E_data[idx(idx_t_left, idx_x_left, idx_y_left)], + p_E_data[idx(idx_t_left, idx_x_left, idx_y_right)], + p_E_data[idx(idx_t_left, idx_x_right, idx_y_left)], + p_E_data[idx(idx_t_left, idx_x_right, idx_y_right)], + p_E_data[idx(idx_t_right, idx_x_left, idx_y_left)], + p_E_data[idx(idx_t_right, idx_x_left, idx_y_right)], + p_E_data[idx(idx_t_right, idx_x_right, idx_y_left)], + p_E_data[idx(idx_t_right, idx_x_right, idx_y_right)], + t, Xp[i], Yp[i])*tmp_e_max; +#endif + } + ); +} + +void +FromTXYEFileLaserProfile::internal_fill_amplitude_nonuniform( + const int idx_t_left, + const int np, + Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, + Real t, Real * AMREX_RESTRICT const amplitude) const +{ + // Copy member variables to tmp copies + // and get pointers to underlying data for GPU. + const auto tmp_e_max = m_common_params.e_max; + const auto tmp_x_min = m_params.x_coords.front(); + const auto tmp_x_max = m_params.x_coords.back(); + const auto tmp_y_min = m_params.y_coords.front(); + const auto tmp_y_max = m_params.y_coords.back(); + const auto p_x_coords = m_params.x_coords.dataPtr(); + const int tmp_x_coords_size = static_cast<int>(m_params.x_coords.size()); + const auto p_y_coords = m_params.y_coords.dataPtr(); + const int tmp_y_coords_size = static_cast<int>(m_params.y_coords.size()); + const auto p_E_data = m_params.E_data.dataPtr(); + const auto tmp_idx_first_time = m_params.first_time_index; + const int idx_t_right = idx_t_left+1; + const auto t_left = m_params.t_coords[idx_t_left]; + const auto t_right = m_params.t_coords[idx_t_right]; + + // Loop through the macroparticle to calculate the proper amplitude + amrex::ParallelFor( + np, + [=] AMREX_GPU_DEVICE (int i) { + //Amplitude is zero if we are out of bounds + if (Xp[i] <= tmp_x_min || Xp[i] >= tmp_x_max){ + amplitude[i] = 0.0_rt; + return; + } +#if (AMREX_SPACEDIM == 3) + if (Yp[i] <= tmp_y_min || Yp[i] >= tmp_y_max){ + amplitude[i] = 0.0_rt; + return; + } +#endif + + //Find indices along x + auto const p_x_right = WarpXUtilAlgo::upper_bound( + p_x_coords, p_x_coords+tmp_x_coords_size, Xp[i]); + const int idx_x_right = p_x_right - p_x_coords; + const int idx_x_left = idx_x_right - 1; + +#if (AMREX_SPACEDIM == 2) + //Interpolate amplitude + const auto idx = [=](int i, int j){ + return (i-tmp_idx_first_time) * tmp_x_coords_size + j; + }; + amplitude[i] = WarpXUtilAlgo::bilinear_interp( + t_left, t_right, + p_x_coords[idx_x_left], p_x_coords[idx_x_right], + p_E_data[idx(idx_t_left, idx_x_left)], + p_E_data[idx(idx_t_left, idx_x_right)], + p_E_data[idx(idx_t_right, idx_x_left)], + p_E_data[idx(idx_t_right, idx_x_right)], + t, Xp[i])*tmp_e_max; + +#elif (AMREX_SPACEDIM == 3) + //Find indices along y + auto const p_y_right = WarpXUtilAlgo::upper_bound( + p_y_coords, p_y_coords+tmp_y_coords_size, Yp[i]); + const int idx_y_right = p_y_right - p_y_coords; + const int idx_y_left = idx_y_right - 1; + + //Interpolate amplitude + const auto idx = [=](int i, int j, int k){ + return + (i-tmp_idx_first_time)*tmp_x_coords_size*tmp_y_coords_size+ + j*tmp_y_coords_size + k; + }; + amplitude[i] = WarpXUtilAlgo::trilinear_interp( + t_left, t_right, + p_x_coords[idx_x_left], p_x_coords[idx_x_right], + p_y_coords[idx_y_left], p_y_coords[idx_y_right], + p_E_data[idx(idx_t_left, idx_x_left, idx_y_left)], + p_E_data[idx(idx_t_left, idx_x_left, idx_y_right)], + p_E_data[idx(idx_t_left, idx_x_right, idx_y_left)], + p_E_data[idx(idx_t_left, idx_x_right, idx_y_right)], + p_E_data[idx(idx_t_right, idx_x_left, idx_y_left)], + p_E_data[idx(idx_t_right, idx_x_left, idx_y_right)], + p_E_data[idx(idx_t_right, idx_x_right, idx_y_left)], + p_E_data[idx(idx_t_right, idx_x_right, idx_y_right)], + t, Xp[i], Yp[i])*tmp_e_max; +#endif + } + ); +} diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp index a0b5dd855..18bdec4a7 100644 --- a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp +++ b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp @@ -72,7 +72,7 @@ GaussianLaserProfile::init ( void GaussianLaserProfile::fill_amplitude ( const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, - Real t, Real * AMREX_RESTRICT const amplitude) + Real t, Real * AMREX_RESTRICT const amplitude) const { Complex I(0,1); // Calculate a few factors which are independent of the macroparticle diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp index 55374c5ea..7fe75cf56 100644 --- a/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp +++ b/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp @@ -35,7 +35,7 @@ HarrisLaserProfile::init ( void HarrisLaserProfile::fill_amplitude ( const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, - Real t, Real * AMREX_RESTRICT const amplitude) + Real t, Real * AMREX_RESTRICT const amplitude) const { // This function uses the Harris function as the temporal profile of the pulse const Real omega0 = diff --git a/Source/Laser/LaserProfilesImpl/Make.package b/Source/Laser/LaserProfilesImpl/Make.package index 32284c4e4..2fef27b9f 100644 --- a/Source/Laser/LaserProfilesImpl/Make.package +++ b/Source/Laser/LaserProfilesImpl/Make.package @@ -1,6 +1,7 @@ CEXE_sources += LaserProfileGaussian.cpp CEXE_sources += LaserProfileHarris.cpp CEXE_sources += LaserProfileFieldFunction.cpp +CEXE_sources += LaserProfileFromTXYEFile.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Laser/LaserProfilesImpl VPATH_LOCATIONS += $(WARPX_HOME)/Source/Laser/LaserProfilesImpl diff --git a/Source/Parser/Make.package b/Source/Parser/Make.package index 5ce02cbda..15115c138 100644 --- a/Source/Parser/Make.package +++ b/Source/Parser/Make.package @@ -5,6 +5,7 @@ CEXE_sources += WarpXParser.cpp CEXE_headers += WarpXParser.H CEXE_headers += GpuParser.H CEXE_sources += GpuParser.cpp +CEXE_headers += WarpXParserWrapper.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parser VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parser diff --git a/Source/Parser/WarpXParserWrapper.H b/Source/Parser/WarpXParserWrapper.H new file mode 100644 index 000000000..2dd7f72c7 --- /dev/null +++ b/Source/Parser/WarpXParserWrapper.H @@ -0,0 +1,35 @@ +#ifndef WARPX_PARSER_WRAPPER_H_ +#define WARPX_PARSER_WRAPPER_H_ + +#include <WarpXParser.H> +#include <AMReX_Gpu.H> +#include <GpuParser.H> +/** + * \brief + * The ParserWrapper struct is constructed to safely use the GpuParser class + * that can essentially be though of as a raw pointer. The GpuParser does + * not have an explicit destructor and the AddPlasma subroutines handle the GpuParser + * in a safe way. The ParserWrapper struct is used to avoid memory leak + * in the EB parser functions. + */ +struct ParserWrapper + : public amrex::Gpu::Managed +{ + ParserWrapper (WarpXParser const& a_parser) noexcept + : m_parser(a_parser) {} + + ~ParserWrapper() { + m_parser.clear(); + } + + AMREX_GPU_HOST_DEVICE + amrex::Real + getField (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + return m_parser(x,y,z); + } + + GpuParser m_parser; +}; + +#endif diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index c577da7f3..e05a64bfe 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -1,4 +1,5 @@ #include <WarpX.H> +#include <WarpXUtil.H> #include <WarpXConst.H> using namespace amrex; @@ -97,10 +98,25 @@ WarpX::MoveWindow (bool move_j) // Shift each component of vector fields (E, B, j) for (int dim = 0; dim < 3; ++dim) { - // Fine grid - shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir, B_external_grid[dim]); - shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir, E_external_grid[dim]); + ParserWrapper *Bfield_parser; + ParserWrapper *Efield_parser; + bool use_Bparser = false; + bool use_Eparser = false; + if (B_ext_grid_s == "parse_b_ext_grid_function") { + use_Bparser = true; + if (dim == 0) Bfield_parser = Bxfield_parser.get(); + if (dim == 1) Bfield_parser = Byfield_parser.get(); + if (dim == 2) Bfield_parser = Bzfield_parser.get(); + } + if (E_ext_grid_s == "parse_e_ext_grid_function") { + use_Eparser = true; + if (dim == 0) Efield_parser = Exfield_parser.get(); + if (dim == 1) Efield_parser = Eyfield_parser.get(); + if (dim == 2) Efield_parser = Ezfield_parser.get(); + } + shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir, B_external_grid[dim], use_Bparser, Bfield_parser); + shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir, E_external_grid[dim], use_Eparser, Efield_parser); if (move_j) { shiftMF(*current_fp[lev][dim], geom[lev], num_shift, dir); } @@ -112,9 +128,9 @@ WarpX::MoveWindow (bool move_j) } if (lev > 0) { - // Coarse grid - shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim]); - shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, E_external_grid[dim]); + // coarse grid + shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim], use_Bparser, Bfield_parser); + shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, E_external_grid[dim], use_Eparser, Efield_parser); shiftMF(*Bfield_aux[lev][dim], geom[lev], num_shift, dir); shiftMF(*Efield_aux[lev][dim], geom[lev], num_shift, dir); if (move_j) { @@ -132,7 +148,7 @@ WarpX::MoveWindow (bool move_j) // Shift scalar component F for dive cleaning if (do_dive_cleaning) { // Fine grid - shiftMF(*F_fp[lev], geom[lev], num_shift, dir); + shiftMF(*F_fp[lev], geom[lev], num_shift, dir); if (do_pml && pml[lev]->ok()) { MultiFab* pml_F = pml[lev]->GetF_fp(); shiftMF(*pml_F, geom[lev], num_shift, dir); @@ -204,7 +220,7 @@ WarpX::MoveWindow (bool move_j) void WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir, - amrex::Real external_field) + amrex::Real external_field, bool useparser, ParserWrapper *field_parser) { BL_PROFILE("WarpX::shiftMF()"); const BoxArray& ba = mf.boxArray(); @@ -246,20 +262,57 @@ WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir, shiftiv[dir] = num_shift; Dim3 shift = shiftiv.dim3(); + const RealBox& real_box = geom.ProbDomain(); + const auto dx = geom.CellSizeArray(); + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif + + for (MFIter mfi(tmpmf); mfi.isValid(); ++mfi ) { auto const& dstfab = mf.array(mfi); auto const& srcfab = tmpmf.array(mfi); const Box& outbox = mfi.fabbox() & adjBox; + if (outbox.ok()) { - AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n, - { - srcfab(i,j,k,n) = external_field; - }); + if (useparser == false) { + AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n, + { + srcfab(i,j,k,n) = external_field; + }); + } else if (useparser == true) { + // index type of the src mf + auto const& mf_IndexType = (tmpmf).ixType(); + IntVect mf_type(AMREX_D_DECL(0,0,0)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + mf_type[idim] = mf_IndexType.nodeCentered(idim); + } + + amrex::ParallelFor (outbox, nc, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + // Compute x,y,z co-ordinates based on index type of mf + Real fac_x = (1.0 - mf_type[0]) * dx[0]*0.5; + Real x = i*dx[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mf_type[1]) * dx[1]*0.5; + Real z = j*dx[1] + real_box.lo(1) + fac_z; +#else + Real fac_y = (1.0 - mf_type[1]) * dx[1]*0.5; + Real y = j*dx[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mf_type[2]) * dx[2]*0.5; + Real z = k*dx[2] + real_box.lo(2) + fac_z; +#endif + srcfab(i,j,k,n) = field_parser->getField(x,y,z); + } + , amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double)*3 + ); + } + } Box dstBox = mf[mfi].box(); diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index 195e309cf..e7b2ef196 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -4,6 +4,8 @@ #include <AMReX_REAL.H> #include <AMReX_Vector.H> #include <AMReX_MultiFab.H> +#include <AMReX_ParmParse.H> +#include <WarpXParser.H> #include <string> @@ -15,14 +17,108 @@ void ConvertLabParamsToBoost(); void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax); +/** + * \brief Parse a string (typically a mathematical expression) from the + * input file and store it into a variable. + * + * \param ParmParse pp used to read the query_string pp.<function>=string + * \param parmparse_string String used to initialize ParmParse + * \param query_string ParmParse.query will look for this string + * \param stored_string variable in which the string to parse is stored + */ +void Store_parserString(amrex::ParmParse &pp, std::string query_string, + std::string& stored_string); + namespace WarpXUtilIO{ - /** - * A helper function to write binary data on disk. - * @param[in] filename where to write - * @param[in] data Vector containing binary data to write on disk - * return true if it succeeds, false otherwise - */ - bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data); +/** + * A helper function to write binary data on disk. + * @param[in] filename where to write + * @param[in] data Vector containing binary data to write on disk + * return true if it succeeds, false otherwise + */ +bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data); +} + +namespace WarpXUtilAlgo{ + +/** \brief Returns a pointer to the first element in the range [first, last) that is greater than val + * + * A re-implementation of the upper_bound algorithm suitable for GPU kernels. + * + * @param first: pointer to left limit of the range to consider + * @param last: pointer to right limit of the range to consider + * @param val: value to compare the elements of [first, last) to + */ +template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE +const T* upper_bound(const T* first, const T* last, const T& val) +{ + const T* it; + size_t count, step; + count = last-first; + while(count>0){ + it = first; + step = count/2; + it += step; + if (!(val<*it)){ + first = ++it; + count -= step + 1; + } + else{ + count = step; + } + } + return first; } +/** \brief Performs a linear interpolation + * + * Performs a linear interpolation at x given the 2 points + * (x0, f0) and (x1, f1) + */ +template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE +T linear_interp(T x0, T x1, T f0, T f1, T x) +{ + return ((x1-x)*f0 + (x-x0)*f1)/(x1-x0); +} + +/** \brief Performs a bilinear interpolation + * + * Performs a bilinear interpolation at (x,y) given the 4 points + * (x0, y0, f00), (x0, y1, f01), (x1, y0, f10), (x1, y1, f11). + */ +template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE +T bilinear_interp(T x0, T x1, T y0, T y1, T f00, T f01, T f10, T f11, T x, T y) +{ + const T fx0 = linear_interp(x0, x1, f00, f10, x); + const T fx1 = linear_interp(x0, x1, f01, f11, x); + return linear_interp(y0, y1, fx0, fx1, y); +} + +/** \brief Performs a trilinear interpolation + * + * Performs a trilinear interpolation at (x,y,z) given the 8 points + * (x0, y0, z0, f000), (x0, y0, z1, f001), (x0, y1, z0, f010), (x0, y1, z1, f011), + * (x1, y0, z0, f100), (x1, y0, z1, f101), (x1, y1, z0, f110), (x1, y1, z1, f111) + */ +template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE +T trilinear_interp(T x0, T x1,T y0, T y1, T z0, T z1, + T f000, T f001, T f010, T f011, T f100, T f101, T f110, T f111, + T x, T y, T z) +{ + const T fxy0 = bilinear_interp( + x0, x1, y0, y1, f000, f010, f100, f110, x, y); + const T fxy1 = bilinear_interp( + x0, x1, y0, y1, f001, f011, f101, f111, x, y); + return linear_interp(z0, z1, fxy0, fxy1, z); +} + +} + +/** +* \brief Initialize a WarpXParser object from a string containing a math expression +* +* \param parse_function String to read to initialize the parser. +*/ +WarpXParser makeParser (std::string const& parse_function); + #endif //WARPX_UTILS_H_ diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 8764a09c6..e9fb958fd 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -164,3 +164,44 @@ namespace WarpXUtilIO{ } } +void Store_parserString(amrex::ParmParse& pp, std::string query_string, + std::string& stored_string) +{ + + char cstr[query_string.size()+1]; + strcpy(cstr, query_string.c_str()); + + std::vector<std::string> f; + pp.getarr(cstr, f); + stored_string.clear(); + for (auto const& s : f) { + stored_string += s; + } + f.clear(); + +} + + +WarpXParser makeParser (std::string const& parse_function) +{ + WarpXParser parser(parse_function); + parser.registerVariables({"x","y","z"}); + ParmParse pp("my_constants"); + std::set<std::string> symbols = parser.symbols(); + symbols.erase("x"); + symbols.erase("y"); + symbols.erase("z"); + for (auto it = symbols.begin(); it != symbols.end(); ) { + Real v; + if (pp.query(it->c_str(), v)) { + parser.setConstant(*it, v); + it = symbols.erase(it); + } else { + ++it; + } + } + for (auto const& s : symbols) { + amrex::Abort("makeParser::Unknown symbol "+s); + } + return parser; +} diff --git a/Source/WarpX.H b/Source/WarpX.H index decde3dc2..3bab73833 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -24,6 +24,7 @@ #include <AMReX_LayoutData.H> #include <AMReX_Interpolater.H> #include <AMReX_FillPatchUtil.H> +#include <WarpXParserWrapper.H> #ifdef _OPENMP # include <omp.h> @@ -73,8 +74,10 @@ public: MultiParticleContainer& GetPartContainer () { return *mypc; } - static void shiftMF(amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, int dir, - amrex::Real external_field = 0.); + static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, + int num_shift, int dir, amrex::Real external_field=0.0, + bool useparser = false, + ParserWrapper *field_parser=nullptr); static void GotoNextLine (std::istream& is); @@ -86,6 +89,28 @@ public: static amrex::Vector<amrex::Real> E_external_grid; static amrex::Vector<amrex::Real> B_external_grid; + // Initialization Type for External E and B + static std::string B_ext_grid_s; + static std::string E_ext_grid_s; + + // Parser for B_external on the grid + static std::string str_Bx_ext_grid_function; + static std::string str_By_ext_grid_function; + static std::string str_Bz_ext_grid_function; + // Parser for E_external on the grid + static std::string str_Ex_ext_grid_function; + static std::string str_Ey_ext_grid_function; + static std::string str_Ez_ext_grid_function; + + // ParserWrapper for B_external on the grid + std::unique_ptr<ParserWrapper> Bxfield_parser; + std::unique_ptr<ParserWrapper> Byfield_parser; + std::unique_ptr<ParserWrapper> Bzfield_parser; + // ParserWrapper for E_external on the grid + std::unique_ptr<ParserWrapper> Exfield_parser; + std::unique_ptr<ParserWrapper> Eyfield_parser; + std::unique_ptr<ParserWrapper> Ezfield_parser; + // Algorithms static long current_deposition_algo; static long charge_deposition_algo; @@ -323,6 +348,49 @@ public: protected: + /** + * \brief + * This function initializes E, B, rho, and F, at all the levels + * of the multifab. rho and F are initialized with 0. + * The E and B fields are initialized using user-defined inputs. + * The initialization type is set using "B_ext_grid_init_style" + * and "E_ext_grid_init_style". The initialization style is set to "default" + * if not explicitly defined by the user, and the E and B fields are + * initialized with E_external_grid and B_external_grid, respectively, each with + * a default value of 0. + * If the initialization type for the E and B field is "constant", + * then, the E and B fields at all the levels are initialized with + * user-defined values for E_external_grid and B_external_grid. + * If the initialization type for B-field is set to + * "parse_B_ext_grid_function", then, the parser is used to read + * Bx_external_grid_function(x,y,z), By_external_grid_function(x,y,z), + * and Bz_external_grid_function(x,y,z). + * Similarly, if the E-field initialization type is set to + * "parse_E_ext_grid_function", then, the parser is used to read + * Ex_external_grid_function(x,y,z), Ey_external_grid_function(x,y,z), + * and Ex_external_grid_function(x,y,z). The parser for the E and B + * initialization assumes that the function has three independent + * variables, at max, namely, x, y, z. However, any number of constants + * can be used in the function used to define the E and B fields on the grid. + */ + void InitLevelData (int lev, amrex::Real time); + + /** + * \brief + * This function initializes the E and B fields on each level + * using the parser and the user-defined function for the external fields. + * The subroutine will parse the x_/y_z_external_grid_function and + * then, the B or E multifab is initialized based on the (x,y,z) position + * on the staggered yee-grid or cell-centered grid. + */ + void InitializeExternalFieldsOnGridUsingParser ( + amrex::MultiFab *mfx, amrex::MultiFab *mfy, amrex::MultiFab *mfz, + ParserWrapper *xfield_parser, ParserWrapper *yfield_parser, + ParserWrapper *zfield_parser, amrex::IntVect x_nodal_flag, + amrex::IntVect y_nodal_flag, amrex::IntVect z_nodal_flag, + const int lev); + + //! Tagging cells for refinement virtual void ErrorEst (int lev, amrex::TagBoxArray& tags, amrex::Real time, int /*ngrow*/) final; @@ -410,7 +478,6 @@ private: void AllocLevelData (int lev, const amrex::BoxArray& new_grids, const amrex::DistributionMapping& new_dmap); - void InitLevelData (int lev, amrex::Real time); void InitFromCheckpoint (); void PostRestart (); @@ -665,4 +732,5 @@ private: int insitu_pin_mesh; }; + #endif diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 377d103d1..48b4bbd55 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -30,6 +30,18 @@ Vector<Real> WarpX::E_external_particle(3, 0.0); Vector<Real> WarpX::E_external_grid(3, 0.0); Vector<Real> WarpX::B_external_grid(3, 0.0); +std::string WarpX::B_ext_grid_s = "default"; +std::string WarpX::E_ext_grid_s = "default"; + +// Parser for B_external on the grid +std::string WarpX::str_Bx_ext_grid_function; +std::string WarpX::str_By_ext_grid_function; +std::string WarpX::str_Bz_ext_grid_function; +// Parser for E_external on the grid +std::string WarpX::str_Ex_ext_grid_function; +std::string WarpX::str_Ey_ext_grid_function; +std::string WarpX::str_Ez_ext_grid_function; + int WarpX::do_moving_window = 0; int WarpX::moving_window_dir = -1; Real WarpX::moving_window_v = std::numeric_limits<amrex::Real>::max(); @@ -309,9 +321,6 @@ WarpX::ReadParameters () pp.queryarr("B_external_particle", B_external_particle); pp.queryarr("E_external_particle", E_external_particle); - pp.queryarr("E_external_grid", E_external_grid); - pp.queryarr("B_external_grid", B_external_grid); - pp.query("do_moving_window", do_moving_window); if (do_moving_window) { |