diff options
36 files changed, 723 insertions, 241 deletions
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index c81188319..202b25d8e 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -24,15 +24,14 @@ First, let us setup your local git repo. Make your own fork of the main on the [WarpX Github page](https://github.com/ECP-WarpX/WarpX), press the fork button. Then, you can execute: ``` -# These 5 first lines are the same as for a standard WarpX install +# These 4 first lines are the same as for a standard WarpX install mkdir warpx_directory cd warpx_directory -git clone https://bitbucket.org/berkeleylab/picsar.git -git clone https://github.com/AMReX-Codes/amrex.git -cd amrex && git checkout development && cd .. # switch to AMReX development branch +git clone --branch master https://bitbucket.org/berkeleylab/picsar.git +git clone --branch development https://github.com/AMReX-Codes/amrex.git # Clone your fork on your local computer. You can get this address on your fork's Github page. -git clone https://github.com/<myGithubUsername>/ECP-WarpX/WarpX.git +git clone --branch dev https://github.com/<myGithubUsername>/ECP-WarpX/WarpX.git cd warpx # Keep track of the main WarpX repo, to remain up-to-date. git remote add upstream https://github.com/ECP-WarpX/WarpX.git diff --git a/Docs/source/building/building.rst b/Docs/source/building/building.rst index fccbb409a..0094f2763 100644 --- a/Docs/source/building/building.rst +++ b/Docs/source/building/building.rst @@ -17,25 +17,9 @@ single directory (e.g. ``warpx_directory``): mkdir warpx_directory cd warpx_directory - git clone https://github.com/ECP-WarpX/WarpX.git - git clone https://bitbucket.org/berkeleylab/picsar.git - git clone https://github.com/AMReX-Codes/amrex.git - -Then switch to the branch ``development`` of AMReX - -:: - - cd amrex/ - git checkout development - cd .. - -and to the branch ``dev`` of WarpX - -:: - - cd WarpX/ - git checkout dev - cd .. + git clone --branch dev https://github.com/ECP-WarpX/WarpX.git + git clone --branch master https://bitbucket.org/berkeleylab/picsar.git + git clone --branch developement https://github.com/AMReX-Codes/amrex.git Basic compilation ----------------- diff --git a/Docs/source/building/summit.rst b/Docs/source/building/summit.rst index 1cda5ad00..44cc5c0ee 100644 --- a/Docs/source/building/summit.rst +++ b/Docs/source/building/summit.rst @@ -11,21 +11,9 @@ correct branch: mkdir warpx_directory cd warpx_directory - git clone https://github.com/ECP-WarpX/WarpX.git - cd WarpX - git checkout dev - cd .. - - git clone https://bitbucket.org/berkeleylab/picsar.git - cd picsar - git checkout master - cd .. - - git clone https://github.com/AMReX-Codes/amrex.git - cd amrex - git checkout development - cd .. - + git clone --branch dev https://github.com/ECP-WarpX/WarpX.git + git clone --branch master https://bitbucket.org/berkeleylab/picsar.git + git clone --branch development https://github.com/AMReX-Codes/amrex.git Then, use the following set of commands to compile: diff --git a/Docs/source/conf.py b/Docs/source/conf.py index 7604708f8..54e533469 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -58,7 +58,7 @@ author = 'WarpX collaboration' # built documents. # # The short X.Y version. -version = '19.04' +version = '19.05' # The full version, including alpha/beta/rc tags. release = '' diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 505d9b31d..32ac12104 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -171,6 +171,10 @@ Particle initialization Whether particle's weight is varied with their radius. This only applies to cylindrical geometry. The only valid value is true. + * ``predefined``: use one of WarpX predefined plasma profiles. It requires additional + arguments ``<species_name>.predefined_profile_name`` and + ``<species_name>.predefined_profile_params`` (see below). + * ``<species_name>.momentum_distribution_type`` (`string`) Distribution of the normalized momentum (`u=p/mc`) for this species. The options are: @@ -206,6 +210,37 @@ Particle initialization * If ``true``, each particle is advanced with the average speed of the species ``vzbar`` until it reaches ``zinject_plane``. +* ``species_name.predefined_profile_name`` (`string`) + Only read of ``<species_name>.electrons.profile`` is `predefined`. + + * If ``parabolic_channel``, the plasma profile is a parabolic profile with linear ramps + at the beginning and the end of the profile. The density is given by + + .. math:: + + n = n_0 n(x,y) n(z) + + with + + .. math:: + + n(x,y) = 1 + 4\frac{x^2+y^2}{k_p^2 R_c^4} + + where :math:`k_p` is the plasma wavenumber associated with density :math:`n_0`. + Here, :math:`n(z)` is a linear up-ramp from :math:`0` to :math:`L_{ramp,up}`, + constant to :math:`1` from :math:`L_{ramp,up}` to :math:`L_{ramp,up} + L_{plateau}` + and a linear down-ramp from :math:`L_{ramp,up} + L_{plateau}` to + :math:`L_{ramp,up} + L_{plateau}+L_{ramp,down}`. All parameters are given + in ``predefined_profile_params``. + +* ``<species_name>.predefined_profile_params`` (list of `float`) + Parameters for the predefined profiles. + + * If ``species_name.predefined_profile_name`` is ``parabolic_channel``, + ``predefined_profile_params`` contains a space-separated list of the + following parameters, in this order: :math:`L_{ramp,up}` :math:`L_{plateau}` + :math:`L_{ramp,down}` :math:`R_c` :math:`n_0` + * ``<species_name>.do_backward_injection`` (`bool`) Inject a backward-propagating beam to reduce the effect of charge-separation fields when running in the boosted frame. See examples. @@ -373,6 +408,25 @@ Laser initialization Temporal chirp at focus. See definition in Akturk et al., Opt Express, vol 12, no 19 (2014). +* ``warpx.num_mirrors`` (`int`) optional (default `0`) + Users can input perfect mirror condition inside the simulation domain. + The number of mirrors is given by ``warpx.num_mirrors``. The mirrors are + orthogonal to the `z` direction. The following parameters are required + when ``warpx.num_mirrors`` is >0. + +* ``warpx.mirror_z`` (list of `float`) required if ``warpx.num_mirrors>0`` + ``z`` location of the front of the mirrors. + +* ``warpx.mirror_z_width`` (list of `float`) required if ``warpx.num_mirrors>0`` + ``z`` width of the mirrors. + +* ``warpx.mirror_z_npoints`` (list of `int`) required if ``warpx.num_mirrors>0`` + In the boosted frame, depending on `gamma_boost`, ``warpx.mirror_z_width`` + can be smaller than the cell size, so that the mirror would not work. This + parameter is the minimum number of points for the mirror. If + ``mirror_z_width < dz/cell_size``, the upper bound of the mirror is increased + so that it contains at least ``mirror_z_npoints``. + Numerics and algorithms ----------------------- diff --git a/Docs/source/visualization/sensei.rst b/Docs/source/visualization/sensei.rst index 59f879e47..aa372b35b 100644 --- a/Docs/source/visualization/sensei.rst +++ b/Docs/source/visualization/sensei.rst @@ -250,13 +250,10 @@ First, log into cori and clone the git repo's. cd $SCRATCH mkdir warpx cd warpx/ - git clone https://github.com/ECP-WarpX/WarpX.git WarpX-libsim - git clone https://github.com/AMReX-Codes/amrex - git clone https://bitbucket.org/berkeleylab/picsar.git - cd amrex/ - git checkout development - cd ../WarpX-libsim - git checkout dev + git clone --branch dev https://github.com/ECP-WarpX/WarpX.git WarpX-libsim + git clone --branch development https://github.com/AMReX-Codes/amrex + git clone --branch master https://bitbucket.org/berkeleylab/picsar.git + cd WarpX-libsim vim GNUmakefile Next, edit the makefile to turn the SENSEI features on. @@ -304,13 +301,10 @@ First, log into cori and clone the git repo's. cd $SCRATCH mkdir warpx cd warpx/ - git clone https://github.com/ECP-WarpX/WarpX.git WarpX-catalyst - git clone https://github.com/AMReX-Codes/amrex - git clone https://bitbucket.org/berkeleylab/picsar.git - cd amrex/ - git checkout development - cd ../WarpX-catalyst - git checkout dev + git clone --branch dev https://github.com/ECP-WarpX/WarpX.git WarpX-catalyst + git clone --branch development https://github.com/AMReX-Codes/amrex + git clone --branch master https://bitbucket.org/berkeleylab/picsar.git + cd WarpX-catalyst vim GNUmakefile Next, edit the makefile to turn the SENSEI features on. diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost index e07e3a87a..3b452851e 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost +++ b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost @@ -48,12 +48,6 @@ warpx.dt_snapshots_lab = 1.6678204759907604e-12 ################################# ############ PLASMA ############# ################################# -my_constants.zmax = 30.e-2 -my_constants.rc = 50.e-6 -my_constants.rmax = 110.e-6 -my_constants.n0 = 3.5e24 -my_constants.kp = 353352. - particles.nspecies = 3 particles.species_names = electrons ions beam particles.use_fdtd_nci_corr = 1 @@ -66,10 +60,12 @@ electrons.num_particles_per_cell_each_dim = 1 1 electrons.momentum_distribution_type = "gaussian" electrons.xmin = -120.e-6 electrons.xmax = 120.e-6 -electrons.zmin = 0.e-6 -electrons.zmax = .003 -electrons.profile = "parse_density_function" -electrons.density_function(x,y,z) = "((x**2+y**2)<rmax**2)*n0*(1+4*(x**2+y**2)/(kp**2*rc**4))" +electrons.zmin = 0.5e-3 +electrons.zmax = .0035 +electrons.profile = "predefined" +electrons.predefined_profile_name = "parabolic_channel" +# predefined_profile_params = z_start ramp_up plateau ramp_down rc n0 +electrons.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 ions.charge = q_e ions.mass = m_p @@ -78,10 +74,12 @@ ions.num_particles_per_cell_each_dim = 1 1 ions.momentum_distribution_type = "gaussian" ions.xmin = -120.e-6 ions.xmax = 120.e-6 -ions.zmin = 0. -ions.zmax = .003 -ions.profile = "parse_density_function" -ions.density_function(x,y,z) = "((x**2+y**2)<rmax**2)*n0*(1+4*(x**2+y**2)/(kp**2*rc**4))" +ions.zmin = 0.5e-3 +ions.zmax = .0035 +ions.profile = "predefined" +ions.predefined_profile_name = "parabolic_channel" +# predefined_profile_params = z_start ramp_up plateau ramp_down rc n0 +ions.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 beam.charge = -q_e beam.mass = m_e @@ -126,5 +124,5 @@ laser1.e_max = 2.e12 # Maximum amplitude of the laser field (in V/m laser1.profile_waist = 45.e-6 # The waist of the laser (in meters) laser1.profile_duration = 20.e-15 # The duration of the laser (in seconds) laser1.profile_t_peak = 40.e-15 # The time at which the laser reaches its peak (in seconds) -laser1.profile_focal_distance = 0.e-6 # Focal distance from the antenna (in meters) +laser1.profile_focal_distance = 0.5e-3 # Focal distance from the antenna (in meters) laser1.wavelength = 0.81e-6 # The wavelength of the laser (in meters) diff --git a/Examples/Tests/Langmuir/inputs.multi.rz.rt b/Examples/Tests/Langmuir/inputs.multi.rz.rt index e5684dd00..6bc84a434 100644 --- a/Examples/Tests/Langmuir/inputs.multi.rz.rt +++ b/Examples/Tests/Langmuir/inputs.multi.rz.rt @@ -11,9 +11,7 @@ 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 = 1 #40 # How often to write plotfiles. "<= 0" means no plotfiles. -warpx.plot_raw_fields = 1 -warpx.plot_raw_fields_guards = 0 +amr.plot_int = 40 # How often to write plotfiles. "<= 0" means no plotfiles. # Geometry geometry.coord_sys = 1 # 0: Cartesian @@ -40,8 +38,11 @@ interpolation.noz = 1 # CFL warpx.cfl = 1.0 +# Having this turned on makes for a more sensitive test +warpx.do_dive_cleaning = 1 + # Parameters for the plasma wave -my_constants.epsilon = 0.001 +my_constants.epsilon = 0.01 my_constants.kp = 266125.0928256017 my_constants.k0 = 314159.2653589793 my_constants.w0 = 5.e-6 diff --git a/Examples/Tests/Langmuir/langmuir_multi_rz_analysis.py b/Examples/Tests/Langmuir/langmuir_multi_rz_analysis.py index 78dc149bf..2ff30eb5e 100755 --- a/Examples/Tests/Langmuir/langmuir_multi_rz_analysis.py +++ b/Examples/Tests/Langmuir/langmuir_multi_rz_analysis.py @@ -18,7 +18,7 @@ from scipy.constants import e, m_e, epsilon_0, c fn = sys.argv[1] # Parameters (these parameters must match the parameters in `inputs.multi.rz.rt`) -epsilon = 0.001 +epsilon = 0.01 n = 2.e24 w0 = 5.e-6 n_osc_z = 2 diff --git a/LICENSE.txt b/LICENSE.txt index 1db7d9ad2..aeeded3a7 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,4 +1,4 @@ -WarpX v19.04 Copyright (c) 2018, The Regents of the University of California, through Lawrence Berkeley National Laboratory, and Lawrence Livermore National Security, LLC, for the operation of Lawrence Livermore National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. +WarpX v19.05 Copyright (c) 2018, The Regents of the University of California, through Lawrence Berkeley National Laboratory, and Lawrence Livermore National Security, LLC, for the operation of Lawrence Livermore National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: diff --git a/Python/setup.py b/Python/setup.py index 7cfd04ace..f4917e7e2 100644 --- a/Python/setup.py +++ b/Python/setup.py @@ -20,7 +20,7 @@ else: package_data = {} setup (name = 'pywarpx', - version = '19.04', + version = '19.05', packages = ['pywarpx'], package_dir = {'pywarpx':'pywarpx'}, description = """Wrapper of WarpX""", diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 60b0b5fa3..e98561be1 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -1,10 +1,10 @@ - #include <cmath> #include <limits> #include <WarpX.H> #include <WarpXConst.H> #include <WarpX_f.H> +#include <WarpXUtil.H> #ifdef WARPX_USE_PY #include <WarpX_py.H> #endif @@ -38,7 +38,7 @@ WarpX::EvolveEM (int numsteps) { Real walltime_beg_step = amrex::second(); - // Start loop on time steps + // Start loop on time steps amrex::Print() << "\nSTEP " << step+1 << " starts ...\n"; #ifdef WARPX_USE_PY if (warpx_py_beforestep) warpx_py_beforestep(); @@ -53,16 +53,16 @@ WarpX::EvolveEM (int numsteps) if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); - // Reset the costs to 0 - for (int lev = 0; lev <= finest_level; ++lev) { - costs[lev]->setVal(0.0); - } + // Reset the costs to 0 + for (int lev = 0; lev <= finest_level; ++lev) { + costs[lev]->setVal(0.0); + } } for (int lev = 0; lev <= finest_level; ++lev) { - // Perform running average of the costs - // (Giving more importance to most recent costs) - (*costs[lev].get()).mult( (1. - 2./load_balance_int) ); + // Perform running average of the costs + // (Giving more importance to most recent costs) + (*costs[lev].get()).mult( (1. - 2./load_balance_int) ); } } @@ -81,8 +81,8 @@ WarpX::EvolveEM (int numsteps) } is_synchronized = false; } else { - // Beyond one step, we have E^{n} and B^{n}. - // Particles have p^{n-1/2} and x^{n}. + // Beyond one step, we have E^{n} and B^{n}. + // Particles have p^{n-1/2} and x^{n}. FillBoundaryE(); FillBoundaryB(); UpdateAuxilaryData(); @@ -97,6 +97,10 @@ WarpX::EvolveEM (int numsteps) amrex::Abort("Unsupported do_subcycling type"); } + if (num_mirrors>0){ + applyMirrors(cur_time); + } + #ifdef WARPX_USE_PY if (warpx_py_beforeEsolve) warpx_py_beforeEsolve(); #endif @@ -105,8 +109,10 @@ WarpX::EvolveEM (int numsteps) UpdateAuxilaryData(); for (int lev = 0; lev <= finest_level; ++lev) { mypc->PushP(lev, 0.5*dt[lev], - *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], - *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); + *Efield_aux[lev][0],*Efield_aux[lev][1], + *Efield_aux[lev][2], + *Bfield_aux[lev][0],*Bfield_aux[lev][1], + *Bfield_aux[lev][2]); } is_synchronized = true; } @@ -118,18 +124,18 @@ WarpX::EvolveEM (int numsteps) ++istep[lev]; } - cur_time += dt[0]; + cur_time += dt[0]; bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && - (insitu_int > 0) && ((step+1) % insitu_int == 0); + (insitu_int > 0) && ((step+1) % insitu_int == 0); bool move_j = is_synchronized || to_make_plot || do_insitu; // 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(move_j); + int num_moved = MoveWindow(move_j); if (max_level == 0) { int num_redistribute_ghost = num_moved + 1; @@ -139,11 +145,11 @@ WarpX::EvolveEM (int numsteps) mypc->Redistribute(); } - bool to_sort = (sort_int > 0) && ((step+1) % sort_int == 0); - if (to_sort) { - amrex::Print() << "re-sorting particles \n"; - mypc->SortParticlesByCell(); - } + bool to_sort = (sort_int > 0) && ((step+1) % sort_int == 0); + if (to_sort) { + amrex::Print() << "re-sorting particles \n"; + mypc->SortParticlesByCell(); + } amrex::Print()<< "STEP " << step+1 << " ends." << " TIME = " << cur_time << " DT = " << dt[0] << "\n"; @@ -153,10 +159,10 @@ WarpX::EvolveEM (int numsteps) << " s; This step = " << walltime_end_step-walltime_beg_step << " s; Avg. per step = " << walltime/(step+1) << " s\n"; - // sync up time - for (int i = 0; i <= max_level; ++i) { - t_new[i] = cur_time; - } + // sync up time + for (int i = 0; i <= max_level; ++i) { + t_new[i] = cur_time; + } if (do_boosted_frame_diagnostic) { std::unique_ptr<MultiFab> cell_centered_data = nullptr; @@ -166,7 +172,7 @@ WarpX::EvolveEM (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - if (to_make_plot || do_insitu) + if (to_make_plot || do_insitu) { FillBoundaryE(); FillBoundaryB(); @@ -178,34 +184,34 @@ WarpX::EvolveEM (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } - last_plot_file_step = step+1; - last_insitu_step = step+1; + last_plot_file_step = step+1; + last_insitu_step = step+1; - if (to_make_plot) - WritePlotFile(); - - if (do_insitu) - UpdateInSitu(); - } + if (to_make_plot) + WritePlotFile(); + if (do_insitu) + UpdateInSitu(); + } - if (check_int > 0 && (step+1) % check_int == 0) { - last_check_file_step = step+1; - WriteCheckPointFile(); - } + if (check_int > 0 && (step+1) % check_int == 0) { + last_check_file_step = step+1; + WriteCheckPointFile(); + } - if (cur_time >= stop_time - 1.e-3*dt[0]) { - max_time_reached = true; - break; - } + if (cur_time >= stop_time - 1.e-3*dt[0]) { + max_time_reached = true; + break; + } #ifdef WARPX_USE_PY if (warpx_py_afterstep) warpx_py_afterstep(); #endif - // End loop on time steps + // End loop on time steps } - bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step && (max_time_reached || istep[0] >= max_step); + bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step + && (max_time_reached || istep[0] >= max_step); bool do_insitu = (insitu_start >= istep[0]) && (insitu_int > 0) && (istep[0] > last_insitu_step) && (max_time_reached || istep[0] >= max_step); @@ -218,8 +224,10 @@ WarpX::EvolveEM (int numsteps) for (int lev = 0; lev <= finest_level; ++lev) { mypc->FieldGather(lev, - *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], - *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); + *Efield_aux[lev][0],*Efield_aux[lev][1], + *Efield_aux[lev][2], + *Bfield_aux[lev][0],*Bfield_aux[lev][1], + *Bfield_aux[lev][2]); } if (write_plot_file) @@ -229,8 +237,9 @@ WarpX::EvolveEM (int numsteps) UpdateInSitu(); } - if (check_int > 0 && istep[0] > last_check_file_step && (max_time_reached || istep[0] >= max_step)) { - WriteCheckPointFile(); + if (check_int > 0 && istep[0] > last_check_file_step && + (max_time_reached || istep[0] >= max_step)) { + WriteCheckPointFile(); } if (do_boosted_frame_diagnostic) { @@ -462,23 +471,23 @@ WarpX::ComputeDt () Real deltat = 0.; if (maxwell_fdtd_solver_id == 0) { - // CFL time step Yee solver + // CFL time step Yee solver #ifdef WARPX_RZ - // Derived semi-analytically by R. Lehe - deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); + // Derived semi-analytically by R. Lehe + deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); #else - deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), - + 1./(dx[1]*dx[1]), - + 1./(dx[2]*dx[2]))) * PhysConst::c ); + deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), + + 1./(dx[1]*dx[1]), + + 1./(dx[2]*dx[2]))) * PhysConst::c ); #endif } else { - // CFL time step CKC solver + // CFL time step CKC solver #if (BL_SPACEDIM == 3) - const Real delta = std::min(dx[0],std::min(dx[1],dx[2])); + const Real delta = std::min(dx[0],std::min(dx[1],dx[2])); #elif (BL_SPACEDIM == 2) - const Real delta = std::min(dx[0],dx[1]); + const Real delta = std::min(dx[0],dx[1]); #endif - deltat = cfl*delta/PhysConst::c; + deltat = cfl*delta/PhysConst::c; } dt.resize(0); dt.resize(max_level+1,deltat); @@ -493,3 +502,61 @@ WarpX::ComputeDt () dt[0] = const_dt; } } + +/* \brief Apply perfect mirror condition inside the box (not at a boundary). + * In practice, set all fields to 0 on a section of the simulation domain + * (as for a perfect conductor with a given thickness). + * The mirror normal direction has to be parallel to the z axis. + */ +void +WarpX::applyMirrors(Real time){ + // Loop over the mirrors + for(int i_mirror=0; i_mirror<num_mirrors; ++i_mirror){ + // Get mirror properties (lower and upper z bounds) + Real z_min = mirror_z[i_mirror]; + Real z_max_tmp = z_min + mirror_z_width[i_mirror]; + // Boost quantities for boosted frame simulations + if (gamma_boost>1){ + z_min = z_min/gamma_boost - PhysConst::c*beta_boost*time; + z_max_tmp = z_max_tmp/gamma_boost - PhysConst::c*beta_boost*time; + } + // Loop over levels + for(int lev=0; lev<=finest_level; lev++){ + // Make sure that the mirror contains at least + // mirror_z_npoints[i_mirror] cells + Real dz = WarpX::CellSize(lev)[2]; + Real z_max = std::max(z_max_tmp, + z_min+mirror_z_npoints[i_mirror]*dz); + // Get fine patch field MultiFabs + MultiFab& Ex = *Efield_fp[lev][0].get(); + MultiFab& Ey = *Efield_fp[lev][1].get(); + MultiFab& Ez = *Efield_fp[lev][2].get(); + MultiFab& Bx = *Bfield_fp[lev][0].get(); + MultiFab& By = *Bfield_fp[lev][1].get(); + MultiFab& Bz = *Bfield_fp[lev][2].get(); + // Set each field to zero between z_min and z_max + NullifyMF(Ex, lev, z_min, z_max); + NullifyMF(Ey, lev, z_min, z_max); + NullifyMF(Ez, lev, z_min, z_max); + NullifyMF(Bx, lev, z_min, z_max); + NullifyMF(By, lev, z_min, z_max); + NullifyMF(Bz, lev, z_min, z_max); + if (lev>0){ + // Get coarse patch field MultiFabs + MultiFab& Ex = *Efield_cp[lev][0].get(); + MultiFab& Ey = *Efield_cp[lev][1].get(); + MultiFab& Ez = *Efield_cp[lev][2].get(); + MultiFab& Bx = *Bfield_cp[lev][0].get(); + MultiFab& By = *Bfield_cp[lev][1].get(); + MultiFab& Bz = *Bfield_cp[lev][2].get(); + // Set each field to zero between z_min and z_max + NullifyMF(Ex, lev, z_min, z_max); + NullifyMF(Ey, lev, z_min, z_max); + NullifyMF(Ez, lev, z_min, z_max); + NullifyMF(Bx, lev, z_min, z_max); + NullifyMF(By, lev, z_min, z_max); + NullifyMF(Bz, lev, z_min, z_max); + } + } + } +} diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index 50127914d..b0ee54095 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,11 +1,12 @@ CEXE_headers += WarpX_Complex.H CEXE_headers += SpectralSolver.H +CEXE_sources += SpectralSolver.cpp CEXE_headers += SpectralFieldData.H CEXE_sources += SpectralFieldData.cpp -CEXE_headers += PsatdAlgorithm.H -CEXE_sources += PsatdAlgorithm.cpp CEXE_headers += SpectralKSpace.H CEXE_sources += SpectralKSpace.cpp +include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H deleted file mode 100644 index acefcc466..000000000 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef WARPX_PSATD_ALGORITHM_H_ -#define WARPX_PSATD_ALGORITHM_H_ - -#include <SpectralKSpace.H> -#include <SpectralFieldData.H> - -/* \brief Class that updates the field in spectral space - * and stores the coefficients of the corresponding update equation. - */ -class PsatdAlgorithm -{ - using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; - - public: - PsatdAlgorithm(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const bool nodal, const amrex::Real dt); - PsatdAlgorithm() = default; // Default constructor - PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; - void pushSpectralFields(SpectralFieldData& f) const; - - private: - // Modified finite-order vectors - KVectorComponent modified_kx_vec, modified_kz_vec; -#if (AMREX_SPACEDIM==3) - KVectorComponent modified_ky_vec; -#endif - SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; -}; - -#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package new file mode 100644 index 000000000..c62c21f44 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += SpectralBaseAlgorithm.H +CEXE_headers += PsatdAlgorithm.H +CEXE_sources += PsatdAlgorithm.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H new file mode 100644 index 000000000..0487e5226 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -0,0 +1,24 @@ +#ifndef WARPX_PSATD_ALGORITHM_H_ +#define WARPX_PSATD_ALGORITHM_H_ + +#include <SpectralBaseAlgorithm.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithm : public SpectralBaseAlgorithm +{ + public: + PsatdAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::Real dt); + // Redefine update equation from base class + virtual void pushSpectralFields(SpectralFieldData& f) const override final; + + private: + SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; +}; + +#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index ada7506c3..37892d35a 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -9,14 +9,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const Real dt) -// Compute and assign the modified k vectors -: modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)), -#if (AMREX_SPACEDIM==3) - modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)), - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal)) -#else - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) -#endif + // Initialize members of base class + : SpectralBaseAlgorithm( spectral_kspace, dm, + norder_x, norder_y, norder_z, nodal ) { const BoxArray& ba = spectral_kspace.spectralspace_ba; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H new file mode 100644 index 000000000..602eb2473 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -0,0 +1,51 @@ +#ifndef WARPX_SPECTRAL_BASE_ALGORITHM_H_ +#define WARPX_SPECTRAL_BASE_ALGORITHM_H_ + +#include <SpectralKSpace.H> +#include <SpectralFieldData.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + * + * `SpectralBaseAlgorithm` is only a base class and cannot be used directly. + * Instead use its subclasses, which implement the specific field update + * equations for a given spectral algorithm. + */ +class SpectralBaseAlgorithm +{ + public: + // Member function that updates the fields in spectral space ; + // meant to be overridden in subclasses + virtual void pushSpectralFields(SpectralFieldData& f) const = 0; + // The destructor should also be a virtual function, so that + // a pointer to subclass of `SpectraBaseAlgorithm` actually + // calls the subclass's destructor. + virtual ~SpectralBaseAlgorithm() {}; + + protected: // Meant to be used in the subclasses + + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; + + // Constructor + SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal) + // Compute and assign the modified k vectors + : modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)), +#if (AMREX_SPACEDIM==3) + modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)), + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal)) +#else + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) +#endif + {}; + + // Modified finite-order vectors + KVectorComponent modified_kx_vec, modified_kz_vec; +#if (AMREX_SPACEDIM==3) + KVectorComponent modified_ky_vec; +#endif +}; + +#endif // WARPX_SPECTRAL_BASE_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 7444452af..d4019a9a3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -1,8 +1,7 @@ #ifndef WARPX_SPECTRAL_SOLVER_H_ #define WARPX_SPECTRAL_SOLVER_H_ -#include <SpectralKSpace.H> -#include <PsatdAlgorithm.H> +#include <SpectralBaseAlgorithm.H> #include <SpectralFieldData.H> /* \brief Top-level class for the electromagnetic spectral solver @@ -14,28 +13,17 @@ class SpectralSolver { public: - // Inline definition of the member functions of `SpectralSolver` + // Inline definition of the member functions of `SpectralSolver`, + // except the constructor (see `SpectralSolver.cpp`) // The body of these functions is short, since the work is done in the // underlying classes `SpectralFieldData` and `PsatdAlgorithm` - /* \brief Initialize the spectral solver */ + // Constructor SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { - // Initialize all structures using the same distribution mapping dm - - // - Initialize k space object (Contains info about the size of - // the spectral space corresponding to each box in `realspace_ba`, - // as well as the value of the corresponding k coordinates) - const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); - // - Initialize the algorithm (coefficients) over this space - algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, - norder_z, nodal, dt ); - // - Initialize arrays for fields in Fourier space + FFT plans - field_data = SpectralFieldData( realspace_ba, k_space, dm ); - }; + const amrex::RealVect dx, const amrex::Real dt ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally @@ -59,14 +47,20 @@ class SpectralSolver /* \brief Update the fields in spectral space, over one timestep */ void pushSpectralFields(){ BL_PROFILE("SpectralSolver::pushSpectralFields"); - algorithm.pushSpectralFields( field_data ); + // Virtual function: the actual function used here depends + // on the sub-class of `SpectralBaseAlgorithm` that was + // initialized in the constructor of `SpectralSolver` + algorithm->pushSpectralFields( field_data ); }; private: SpectralFieldData field_data; // Store field in spectral space // and perform the Fourier transforms - PsatdAlgorithm algorithm; // Contains Psatd coefficients - // and field update equation + std::unique_ptr<SpectralBaseAlgorithm> algorithm; + // Defines field update equation in spectral space, + // and the associated coefficients. + // SpectralBaseAlgorithm is a base class ; this pointer is meant + // to point an instance of a *sub-class* defining a specific algorithm }; #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp new file mode 100644 index 000000000..c21c3cfb1 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -0,0 +1,35 @@ +#include <SpectralKSpace.H> +#include <SpectralSolver.H> +#include <PsatdAlgorithm.H> + +/* \brief Initialize the spectral Maxwell solver + * + * This function selects the spectral algorithm to be used, allocates the + * corresponding coefficients for the discretized field update equation, + * and prepares the structures that store the fields in spectral space. + */ +SpectralSolver::SpectralSolver( + const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::RealVect dx, const amrex::Real dt ) { + + // Initialize all structures using the same distribution mapping dm + + // - Initialize k space object (Contains info about the size of + // the spectral space corresponding to each box in `realspace_ba`, + // as well as the value of the corresponding k coordinates) + const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); + + // - Select the algorithm depending on the input parameters + // Initialize the corresponding coefficients over k space + // TODO: Add more algorithms + selection depending on input parameters + // For the moment, this only uses the standard PsatdAlgorithm + algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + + // - Initialize arrays for fields in spectral space + FFT plans + field_data = SpectralFieldData( realspace_ba, k_space, dm ); + +}; diff --git a/Source/FortranInterface/WarpX_f.F90 b/Source/FortranInterface/WarpX_f.F90 index ccd6dd48a..d7f7a8053 100644 --- a/Source/FortranInterface/WarpX_f.F90 +++ b/Source/FortranInterface/WarpX_f.F90 @@ -186,6 +186,41 @@ contains end subroutine warpx_compute_dive_2d + subroutine warpx_compute_dive_rz (lo, hi, dive, dlo, dhi, & + Ex, xlo, xhi, Ey, ylo, yhi, Ez, zlo, zhi, dx, rmin) & + bind(c, name='warpx_compute_dive_rz') + integer, intent(in) :: lo(2),hi(2),dlo(2),dhi(2),xlo(2),xhi(2),ylo(2),yhi(2),zlo(2),zhi(2) + real(amrex_real), intent(in) :: dx(3), rmin + real(amrex_real), intent(in ) :: Ex (xlo(1):xhi(1),xlo(2):xhi(2)) + real(amrex_real), intent(in ) :: Ey (ylo(1):yhi(1),ylo(2):yhi(2)) + real(amrex_real), intent(in ) :: Ez (zlo(1):zhi(1),zlo(2):zhi(2)) + real(amrex_real), intent(inout) :: dive(dlo(1):dhi(1),dlo(2):dhi(2)) + + integer :: i,k + real(amrex_real) :: dxinv(3) + real(amrex_real) :: ru, rd + + dxinv = 1.d0/dx + + do k = lo(2), hi(2) + do i = lo(1), hi(1) + if (i == 0 .and. rmin == 0.) then + ! the bulk equation diverges on axis + ! (due to the 1/r terms). The following expressions regularize + ! these divergences. + dive(i,k) = 4.*dxinv(1) * Ex(i,k) & + + dxinv(3) * (Ez(i,k) - Ez(i,k-1)) + else + ru = 1.d0 + 0.5d0/(rmin*dxinv(1) + i) + rd = 1.d0 - 0.5d0/(rmin*dxinv(1) + i) + dive(i,k) = dxinv(1) * (ru*Ex(i,k) - rd*Ex(i-1,k)) & + + dxinv(3) * (Ez(i,k) - Ez(i,k-1)) + end if + end do + end do + end subroutine warpx_compute_dive_rz + + subroutine warpx_sync_current_2d (lo, hi, crse, clo, chi, fine, flo, fhi, dir) & bind(c, name='warpx_sync_current_2d') integer, intent(in) :: lo(2), hi(2), flo(2), fhi(2), clo(2), chi(2), dir diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index f246f9f54..3d92b7651 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -48,7 +48,6 @@ #elif (AMREX_SPACEDIM == 2) #define WRPX_COMPUTE_DIVB warpx_compute_divb_2d -#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d #define WRPX_SYNC_CURRENT warpx_sync_current_2d #define WRPX_SYNC_RHO warpx_sync_rho_2d @@ -74,6 +73,12 @@ #define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs #define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_2d +#ifdef WARPX_RZ +#define WRPX_COMPUTE_DIVE warpx_compute_dive_rz +#else +#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d +#endif + #endif #ifdef __cplusplus @@ -102,6 +107,12 @@ extern "C" const long* nox, const long* noy,const long* noz, const long* lvect, const long* charge_depo_algo); + // Charge deposition finalize for RZ + void warpx_charge_deposition_rz_volume_scaling( + amrex::Real* rho, const long* rho_ng, const int* rho_ntot, + const amrex::Real* rmin, + const amrex::Real* dr); + // Current deposition void warpx_current_deposition( amrex::Real* jx, const long* jx_ng, const int* jx_ntot, @@ -362,7 +373,11 @@ extern "C" const BL_FORT_FAB_ARG_ANYD(ex), const BL_FORT_FAB_ARG_ANYD(ey), const BL_FORT_FAB_ARG_ANYD(ez), - const amrex::Real* dx); + const amrex::Real* dx +#ifdef WARPX_RZ + ,const amrex::Real* rmin +#endif + ); void WRPX_PUSH_PML_BVEC(const int* xlo, const int* xhi, const int* ylo, const int* yhi, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 151342ff5..6d6a2e411 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -18,7 +18,8 @@ #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2drz_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz -#define WRPX_PXR_RZ_VOLUME_SCALING apply_rz_volume_scaling +#define WRPX_PXR_RZ_VOLUME_SCALING_RHO apply_rz_volume_scaling_rho +#define WRPX_PXR_RZ_VOLUME_SCALING_J apply_rz_volume_scaling_j #define WRPX_PXR_PUSH_BVEC pxrpush_emrz_bvec #define WRPX_PXR_PUSH_EVEC pxrpush_emrz_evec @@ -227,9 +228,15 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! Dimension 2 #elif (AMREX_SPACEDIM==2) +#ifdef WARPX_RZ + logical(pxr_logical) :: l_2drz = .TRUE._c_long +#else + logical(pxr_logical) :: l_2drz = .FALSE._c_long +#endif + CALL pxr_depose_rho_n_2dxz(rho,np,xp,yp,zp,w,q,xmin,zmin,dx,dz,nx,nz,& nxguard,nzguard,nox,noz, & - .TRUE._c_long, .FALSE._c_long, .FALSE._c_long, 0_c_long) + .TRUE._c_long, .FALSE._c_long, l_2drz, 0_c_long) #endif @@ -238,6 +245,44 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! _________________________________________________________________ !> !> @brief + !> Applies the inverse volume scaling for RZ charge deposition + !> + !> @details + !> The scaling is done for both single mode (FDTD) and + !> multi mode (spectral) (todo) + ! + !> @param[inout] rho charge array + !> @param[in] rmin tile grid minimum radius + !> @param[in] dr radial space discretization steps + !> @param[in] nx,ny,nz number of cells + !> @param[in] nxguard,nyguard,nzguard number of guard cells + !> + subroutine warpx_charge_deposition_rz_volume_scaling(rho,rho_ng,rho_ntot,rmin,dr) & + bind(C, name="warpx_charge_deposition_rz_volume_scaling") + + integer, intent(in) :: rho_ntot(AMREX_SPACEDIM) + integer(c_long), intent(in) :: rho_ng + real(amrex_real), intent(IN OUT):: rho(*) + real(amrex_real), intent(IN) :: rmin, dr + + integer(c_long) :: type_rz_depose = 1 + + ! Compute the number of valid cells and guard cells + integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) + rho_nvalid = rho_ntot - 2*rho_ng + rho_nguards = rho_ng + +#ifdef WARPX_RZ + CALL WRPX_PXR_RZ_VOLUME_SCALING_RHO( & + rho,rho_nguards,rho_nvalid, & + rmin,dr,type_rz_depose) +#endif + + end subroutine warpx_charge_deposition_rz_volume_scaling + + ! _________________________________________________________________ + !> + !> @brief !> Main subroutine for the current deposition !> !> @details @@ -353,7 +398,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jz_nguards = jz_ng #ifdef WARPX_RZ - CALL WRPX_PXR_RZ_VOLUME_SCALING( & + CALL WRPX_PXR_RZ_VOLUME_SCALING_J( & jx,jx_nguards,jx_nvalid, & jy,jy_nguards,jy_nvalid, & jz,jz_nguards,jz_nvalid, & @@ -413,7 +458,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n END SELECT !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) #elif (AMREX_SPACEDIM == 2) CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) @@ -498,7 +543,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) #elif (AMREX_SPACEDIM == 2) CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index b825924c6..edcf402c9 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -1,4 +1,5 @@ CEXE_sources += CustomDensityProb.cpp +CEXE_sources += PlasmaProfiles.cpp CEXE_sources += WarpXInitData.cpp CEXE_sources += CustomMomentumProb.cpp CEXE_sources += PlasmaInjector.cpp diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index f6b76e8b6..cd5802a55 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -10,6 +10,8 @@ #include "AMReX_ParmParse.H" #include "AMReX_Utility.H" +enum class predefined_profile_flag { Null, parabolic_channel }; + /// /// PlasmaDensityProfile describes how the charge density /// is set in particle initialization. Subclasses must define a @@ -59,6 +61,24 @@ private: }; /// +/// This describes predefined density distributions. +/// +class PredefinedDensityProfile : public PlasmaDensityProfile +{ +public: + PredefinedDensityProfile(const std::string& species_name); + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const override; + amrex::Real ParabolicChannel(amrex::Real x, + amrex::Real y, + amrex::Real z) const; +private: + predefined_profile_flag which_profile = predefined_profile_flag::Null; + amrex::Vector<amrex::Real> params; +}; + +/// /// This describes a density function parsed in the input file. /// class ParseDensityProfile : public PlasmaDensityProfile diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 0da9318de..52b5d8b61 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -70,6 +70,17 @@ CustomDensityProfile::CustomDensityProfile(const std::string& species_name) pp.getarr("custom_profile_params", params); } +PredefinedDensityProfile::PredefinedDensityProfile(const std::string& species_name) +{ + ParmParse pp(species_name); + std::string which_profile_s; + pp.getarr("predefined_profile_params", params); + pp.query("predefined_profile_name", which_profile_s); + if (which_profile_s == "parabolic_channel"){ + which_profile = predefined_profile_flag::parabolic_channel; + } +} + ParseDensityProfile::ParseDensityProfile(std::string parse_density_function) : _parse_density_function(parse_density_function) { @@ -333,6 +344,8 @@ void PlasmaInjector::parseDensity(ParmParse pp){ rho_prof.reset(new ConstantDensityProfile(density)); } else if (rho_prof_s == "custom") { rho_prof.reset(new CustomDensityProfile(species_name)); + } else if (rho_prof_s == "predefined") { + rho_prof.reset(new PredefinedDensityProfile(species_name)); } else if (rho_prof_s == "parse_density_function") { pp.get("density_function(x,y,z)", str_density_function); rho_prof.reset(new ParseDensityProfile(str_density_function)); diff --git a/Source/Initialization/PlasmaProfiles.cpp b/Source/Initialization/PlasmaProfiles.cpp new file mode 100644 index 000000000..d9d207f7e --- /dev/null +++ b/Source/Initialization/PlasmaProfiles.cpp @@ -0,0 +1,41 @@ +#include <PlasmaInjector.H> +#include <cmath> +#include <iostream> +#include <WarpXConst.H> + +using namespace amrex; + +Real PredefinedDensityProfile::getDensity(Real x, Real y, Real z) const { + Real n; + if ( which_profile == predefined_profile_flag::parabolic_channel ) { + n = ParabolicChannel(x,y,z); + } + return n; +} + +/// +/// plateau between linear upramp and downramp, and parab transverse profile +/// +Real PredefinedDensityProfile::ParabolicChannel(Real x, Real y, Real z) const { + // params = [z_start ramp_up plateau ramp_down rc n0] + Real z_start = params[0]; + Real ramp_up = params[1]; + Real plateau = params[2]; + Real ramp_down = params[3]; + Real rc = params[4]; + Real n0 = params[5]; + Real n; + Real kp = PhysConst::q_e/PhysConst::c*sqrt( n0/(PhysConst::m_e*PhysConst::ep0) ); + + if ((z-z_start)>=0 and (z-z_start)<ramp_up ) { + n = (z-z_start)/ramp_up; + } else if ((z-z_start)>=ramp_up and (z-z_start)<ramp_up+plateau ) { + n = 1; + } else if ((z-z_start)>=ramp_up+plateau and (z-z_start)<ramp_up+plateau+ramp_down) { + n = 1-((z-z_start)-ramp_up-plateau)/ramp_down; + } else { + n = 0; + } + n *= n0*(1+4*(x*x+y*y)/(kp*kp*std::pow(rc,4))); + return n; +} diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package index cbb1b5234..3d1fcf1da 100644 --- a/Source/Parallelization/Make.package +++ b/Source/Parallelization/Make.package @@ -1,5 +1,6 @@ CEXE_sources += WarpXComm.cpp CEXE_sources += WarpXRegrid.cpp +CEXE_headers += WarpXSumGuardCells.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parallelization diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 28971eb0c..348b31c8b 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1,5 +1,6 @@ #include <WarpX.H> #include <WarpX_f.H> +#include <WarpXSumGuardCells.H> #include <AMReX_FillPatchUtil_F.H> @@ -355,7 +356,8 @@ WarpX::SyncCurrent () { BL_PROFILE("SyncCurrent()"); - // Restrict fine patch current onto the coarse patch, before fine patch SumBoundary + // Restrict fine patch current onto the coarse patch, before + // summing the guard cells of the fine patch for (int lev = 1; lev <= finest_level; ++lev) { current_cp[lev][0]->setVal(0.0); @@ -391,8 +393,9 @@ WarpX::SyncCurrent () bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); // Then swap j_fp and current_fp std::swap(j_fp[lev][idim], current_fp[lev][idim]); - // At this point, current_fp may have false values close to the - // edges of each FAB. This will be solved with a SumBoundary later. + // At this point, current_fp may have false values close to the + // edges of each FAB. This will be solved with when summing + // the guard cells later. // j_fp contains the exact MultiFab current_fp before this step. } } @@ -427,11 +430,11 @@ WarpX::SyncCurrent () { const auto& period = Geom(lev).periodicity(); // When using a bilinear filter with many passes, current_fp has - // temporarily more ghost cells here, so that its value inside + // temporarily more ghost cells here, so that its value inside // the domain is correct at the end of this stage. - current_fp[lev][0]->SumBoundary(period); - current_fp[lev][1]->SumBoundary(period); - current_fp[lev][2]->SumBoundary(period); + WarpXSumGuardCells(*(current_fp[lev][0]),period); + WarpXSumGuardCells(*(current_fp[lev][1]),period); + WarpXSumGuardCells(*(current_fp[lev][2]),period); } // Add fine level's coarse patch to coarse level's fine patch @@ -461,9 +464,9 @@ WarpX::SyncCurrent () for (int lev = 1; lev <= finest_level; ++lev) { const auto& cperiod = Geom(lev-1).periodicity(); - current_cp[lev][0]->SumBoundary(cperiod); - current_cp[lev][1]->SumBoundary(cperiod); - current_cp[lev][2]->SumBoundary(cperiod); + WarpXSumGuardCells(*(current_cp[lev][0]),cperiod); + WarpXSumGuardCells(*(current_cp[lev][1]),cperiod); + WarpXSumGuardCells(*(current_cp[lev][2]),cperiod); } if (WarpX::use_filter) { @@ -476,7 +479,7 @@ WarpX::SyncCurrent () std::swap(j_fp[lev][idim], current_fp[lev][idim]); // Then copy the interior of j_fp to current_fp. MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, 1, 0); - // current_fp has right number of ghost cells and + // current_fp has right number of ghost cells and // correct filtered values here. } } @@ -556,7 +559,8 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, { if (!rhof[0]) return; - // Restrict fine patch onto the coarse patch, before fine patch SumBoundary + // Restrict fine patch onto the coarse patch, + // before summing the guard cells of the fine patch for (int lev = 1; lev <= finest_level; ++lev) { rhoc[lev]->setVal(0.0); @@ -607,7 +611,7 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, for (int lev = 0; lev <= finest_level; ++lev) { const auto& period = Geom(lev).periodicity(); - rhof[lev]->SumBoundary(period); + WarpXSumGuardCells( *(rhof[lev]), period, 0, rhof[lev]->nComp() ); } // Add fine level's coarse patch to coarse level's fine patch @@ -631,7 +635,7 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { const auto& cperiod = Geom(lev-1).periodicity(); - rhoc[lev]->SumBoundary(cperiod); + WarpXSumGuardCells( *(rhoc[lev]), cperiod, 0, rhoc[lev]->nComp() ); } if (WarpX::use_filter) { @@ -729,10 +733,9 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *j[idim]); - jf.SumBoundary(period); - MultiFab::Copy(*j[idim], jf, 0, 0, 1, 0); + WarpXSumGuardCells(*(j[idim]), jf, period); } else { - j[idim]->SumBoundary(period); + WarpXSumGuardCells(*(j[idim]), period); } } } @@ -780,8 +783,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) MultiFab::Add(jfb, jfc, 0, 0, 1, ng); mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); - jfc.SumBoundary(period); - MultiFab::Copy(*current_cp[lev+1][idim], jfc, 0, 0, 1, 0); + WarpXSumGuardCells(*current_cp[lev+1][idim], jfc, period); } else if (use_filter) // but no buffer { @@ -792,8 +794,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_cp[lev+1][idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]); mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); - jf.SumBoundary(period); - MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0); + WarpXSumGuardCells(*current_cp[lev+1][idim], jf, period); } else if (current_buf[lev+1][idim]) // but no filter { @@ -803,14 +804,14 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, 1, current_buf[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); - current_cp[lev+1][idim]->SumBoundary(period); + WarpXSumGuardCells(*(current_cp[lev+1][idim]), period); } else // no filter, no buffer { mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); - current_cp[lev+1][idim]->SumBoundary(period); + WarpXSumGuardCells(*(current_cp[lev+1][idim]), period); } MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } @@ -840,10 +841,9 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp); - rf.SumBoundary(period); - MultiFab::Copy(*r, rf, 0, icomp, ncomp, 0); + WarpXSumGuardCells(*r, rf, period, icomp, ncomp ); } else { - r->SumBoundary(icomp, ncomp, period); + WarpXSumGuardCells(*r, period, icomp, ncomp); } } @@ -885,9 +885,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) MultiFab::Add(rhofb, rhofc, 0, 0, ncomp, ng); mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); - - rhofc.SumBoundary(period); - MultiFab::Copy(*rho_cp[lev+1], rhofc, 0, 0, ncomp, 0); + WarpXSumGuardCells( *rho_cp[lev+1], rhofc, period, icomp, ncomp ); } else if (use_filter) // but no buffer { @@ -896,8 +894,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp); mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); - rf.SumBoundary(0, ncomp, period); - MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); + WarpXSumGuardCells( *rho_cp[lev+1], rf, period, icomp, ncomp ); } else if (charge_buf[lev+1]) // but no filter { @@ -908,14 +905,14 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) ncomp, charge_buf[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); - rho_cp[lev+1]->SumBoundary(icomp, ncomp, period); + WarpXSumGuardCells(*(rho_cp[lev+1]), period, icomp, ncomp); } else // no filter, no buffer { mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); - rho_cp[lev+1]->SumBoundary(icomp, ncomp, period); + WarpXSumGuardCells(*(rho_cp[lev+1]), period, icomp, ncomp); } ApplyFilterandSumBoundaryRho(lev, PatchType::fine, icomp, ncomp); MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0); diff --git a/Source/Parallelization/WarpXSumGuardCells.H b/Source/Parallelization/WarpXSumGuardCells.H new file mode 100644 index 000000000..24ad1b80f --- /dev/null +++ b/Source/Parallelization/WarpXSumGuardCells.H @@ -0,0 +1,61 @@ +#ifndef WARPX_SUM_GUARD_CELLS_H_ +#define WARPX_SUM_GUARD_CELLS_H_ + +#include <AMReX_MultiFab.H> + +/* \brief Sum the values of `mf`, where the different boxes overlap + * (i.e. in the guard cells) + * + * This is typically called for the sources of the Maxwell equations (J/rho) + * after deposition from the macroparticles. + * + * - When WarpX is compiled with a finite-difference scheme: this only + * updates the *valid* cells of `mf` + * - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this + * updates both the *valid* cells and *guard* cells. (This is because a + * spectral solver requires the value of the sources over a large stencil.) + */ +void +WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period, + const int icomp=0, const int ncomp=1){ +#ifdef WARPX_USE_PSATD + // Update both valid cells and guard cells + const amrex::IntVect n_updated_guards = mf.nGrowVect(); +#else + // Update only the valid cells + const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector(); +#endif + mf.SumBoundary(icomp, ncomp, n_updated_guards, period); +} + +/* \brief Sum the values of `src` where the different boxes overlap + * (i.e. in the guard cells) and copy them into `dst` + * + * This is typically called for the sources of the Maxwell equations (J/rho) + * after deposition from the macroparticles + filtering. + * + * - When WarpX is compiled with a finite-difference scheme: this only + * updates the *valid* cells of `dst` + * - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this + * updates both the *valid* cells and *guard* cells. (This is because a + * spectral solver requires the value of the sources over a large stencil.) + * + * Note: `i_comp` is the component where the results will be stored in `dst`; + * The component from which we copy in `src` is always 0. + */ +void +WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src, + const amrex::Periodicity& period, + const int icomp=0, const int ncomp=1){ +#ifdef WARPX_USE_PSATD + // Update both valid cells and guard cells + const amrex::IntVect n_updated_guards = dst.nGrowVect(); +#else + // Update only the valid cells + const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector(); +#endif + src.SumBoundary(icomp, ncomp, n_updated_guards, period); + amrex::Copy( dst, src, 0, icomp, ncomp, n_updated_guards ); +} + +#endif // WARPX_SUM_GUARD_CELLS_H_ diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 695faaa62..80f3882a0 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -641,6 +641,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); #ifndef AMREX_USE_GPU @@ -696,6 +701,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &cxyzmin_tile[0], &cdx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); #ifndef AMREX_USE_GPU @@ -852,6 +862,12 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + long ngRho = WarpX::nox; + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif #ifdef _OPENMP rhofab.atomicAdd(local_rho); diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index a973c2283..39fded8e8 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -1,7 +1,11 @@ #include <AMReX_REAL.H> #include <AMReX_Vector.H> +#include <AMReX_MultiFab.H> void ReadBoostedFrameParameters(amrex::Real& gamma_boost, amrex::Real& beta_boost, amrex::Vector<int>& boost_direction); void ConvertLabParamsToBoost(); + +void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, + amrex::Real zmax); diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4a884330a..4ec7ebb51 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -3,6 +3,7 @@ #include <WarpXUtil.H> #include <WarpXConst.H> #include <AMReX_ParmParse.H> +#include <WarpX.H> using namespace amrex; @@ -95,3 +96,44 @@ void ConvertLabParamsToBoost() pp_wpx.addarr("fine_tag_hi", fine_tag_hi); } } + +/* \brief Function that sets the value of MultiFab MF to zero for z between + * zmin and zmax. + */ +void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax){ + BL_PROFILE("WarpX::NullifyMF()"); +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for(amrex::MFIter mfi(mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi){ + const amrex::Box& bx = mfi.tilebox(); + // Get box lower and upper physical z bound, and dz + const amrex::Real zmin_box = WarpX::LowerCorner(bx, lev)[2]; + const amrex::Real zmax_box = WarpX::UpperCorner(bx, lev)[2]; + amrex::Real dz = WarpX::CellSize(lev)[2]; + // Get box lower index in the z direction +#if (AMREX_SPACEDIM==3) + const int lo_ind = bx.loVect()[2]; +#else + const int lo_ind = bx.loVect()[1]; +#endif + // Check if box intersect with [zmin, zmax] + if ( (zmin>zmin_box && zmin<=zmax_box) || + (zmax>zmin_box && zmax<=zmax_box) ){ + Array4<Real> arr = mf[mfi].array(); + // Set field to 0 between zmin and zmax + ParallelFor(bx, + [=] 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 + const Real z_gridpoint = zmin_box+(j-lo_ind)*dz; +#endif + if ( (z_gridpoint >= zmin) && (z_gridpoint < zmax) ) { + arr(i,j,k) = 0.; + }; + } + ); + } + } +} diff --git a/Source/WarpX.H b/Source/WarpX.H index 583820646..6cec8e5be 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -144,6 +144,13 @@ public: static amrex::IntVect filter_npass_each_dir; BilinearFilter bilinear_filter; + static int num_mirrors; + amrex::Vector<amrex::Real> mirror_z; + amrex::Vector<amrex::Real> mirror_z_width; + amrex::Vector<int> mirror_z_npoints; + + void applyMirrors(amrex::Real time); + void ComputeDt (); int MoveWindow (bool move_j); void UpdatePlasmaInjectionPosition (amrex::Real dt); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 532858556..47ead98df 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -50,6 +50,8 @@ bool WarpX::use_filter = false; bool WarpX::serialize_ics = false; bool WarpX::refine_plasma = false; +int WarpX::num_mirrors = 0; + int WarpX::sort_int = -1; bool WarpX::do_boosted_frame_diagnostic = false; @@ -355,6 +357,16 @@ WarpX::ReadParameters () filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; #endif + pp.query("num_mirrors", num_mirrors); + if (num_mirrors>0){ + mirror_z.resize(num_mirrors); + pp.getarr("mirror_z", mirror_z, 0, num_mirrors); + mirror_z_width.resize(num_mirrors); + pp.getarr("mirror_z_width", mirror_z_width, 0, num_mirrors); + mirror_z_npoints.resize(num_mirrors); + pp.getarr("mirror_z_npoints", mirror_z_npoints, 0, num_mirrors); + } + pp.query("serialize_ics", serialize_ics); pp.query("refine_plasma", refine_plasma); pp.query("do_dive_cleaning", do_dive_cleaning); @@ -975,12 +987,19 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, for (MFIter mfi(divE, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); +#ifdef WARPX_RZ + const Real xmin = bx.smallEnd(0)*dx[0]; +#endif WRPX_COMPUTE_DIVE(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_N_ANYD(divE[mfi],dcomp), BL_TO_FORTRAN_ANYD((*E[0])[mfi]), BL_TO_FORTRAN_ANYD((*E[1])[mfi]), BL_TO_FORTRAN_ANYD((*E[2])[mfi]), - dx.data()); + dx.data() +#ifdef WARPX_RZ + ,&xmin +#endif + ); } } @@ -995,12 +1014,19 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, for (MFIter mfi(divE, true); mfi.isValid(); ++mfi) { Box bx = mfi.growntilebox(ngrow); +#ifdef WARPX_RZ + const Real xmin = bx.smallEnd(0)*dx[0]; +#endif WRPX_COMPUTE_DIVE(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_N_ANYD(divE[mfi],dcomp), BL_TO_FORTRAN_ANYD((*E[0])[mfi]), BL_TO_FORTRAN_ANYD((*E[1])[mfi]), BL_TO_FORTRAN_ANYD((*E[2])[mfi]), - dx.data()); + dx.data() +#ifdef WARPX_RZ + ,&xmin +#endif + ); } } diff --git a/run_test.sh b/run_test.sh index c0cececb4..a2fec5cd6 100755 --- a/run_test.sh +++ b/run_test.sh @@ -22,9 +22,8 @@ mkdir warpx cp -r ../* warpx # Clone PICSAR and AMReX -git clone https://github.com/AMReX-Codes/amrex.git -cd amrex ; git checkout development ; cd .. -git clone https://bitbucket.org/berkeleylab/picsar.git +git clone --branch development https://github.com/AMReX-Codes/amrex.git +git clone --branch master https://bitbucket.org/berkeleylab/picsar.git # Clone the AMReX regression test utility git clone https://github.com/RemiLehe/regression_testing.git |