diff options
99 files changed, 2459 insertions, 2629 deletions
diff --git a/Docs/source/building/building.rst b/Docs/source/building/building.rst index 33357a4b6..c9c053c95 100644 --- a/Docs/source/building/building.rst +++ b/Docs/source/building/building.rst @@ -106,5 +106,5 @@ Advanced building instructions spectral cori - summitdev + summit spack diff --git a/Docs/source/building/summit.rst b/Docs/source/building/summit.rst new file mode 100644 index 000000000..1cda5ad00 --- /dev/null +++ b/Docs/source/building/summit.rst @@ -0,0 +1,61 @@ +Building WarpX for Summit (OLCF) +================================ + +For the `Summit cluster +<https://www.olcf.ornl.gov/summit/>`__ at OLCF, +use the following commands to download the source code, and switch to the +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 .. + + +Then, use the following set of commands to compile: + +:: + + module load pgi + module load cuda + make -j 4 USE_GPU=TRUE COMP=pgi + +In order to submit a simulation, create a file `submission_script` with +the following text (replace bracketed variables): + +:: + + #!/bin/bash + #BSUB -J <jobName> + #BSUB -W <requestedTime> + #BSUB -nnodes <numberOfNodes> + #BSUB -P <accountNumber> + + module load pgi + module load cuda + + omp=1 + export OMP_NUM_THREADS=${omp} + jsrun -n <numberOfNodes> -a 6 -g 6 -c 6 --bind=packed:${omp} --smpiargs="-gpu" <warpxExecutable> <inputScript> + + +Then run + +:: + + bsub submission_script diff --git a/Docs/source/building/summitdev.rst b/Docs/source/building/summitdev.rst deleted file mode 100644 index c5a3a7034..000000000 --- a/Docs/source/building/summitdev.rst +++ /dev/null @@ -1,36 +0,0 @@ -Building WarpX for Summit-dev (OLCF) -==================================== - -For the `Summit-dev cluster -<https://www.olcf.ornl.gov/tag/summitdev/>`__ at OLCF, -use the following commands to download the source code, and switch to the -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 gpu - cd .. - - git clone https://github.com/AMReX-Codes/amrex.git - cd amrex - git checkout development - cd .. - - -Then, use the following set of commands to compile: - -:: - - module load pgi - module load cuda - make -j 4 USE_GPU=TRUE COMP=pgi diff --git a/Docs/source/conf.py b/Docs/source/conf.py index ce6f436d3..d5165bd17 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 = '0.1.0' +version = '19.02' # 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 e0a0ac2b2..96ec75c05 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -208,6 +208,15 @@ Particle initialization Inject a backward-propagating beam to reduce the effect of charge-separation fields when running in the boosted frame. See examples. +* ``<species_name>.do_splitting`` (`bool`) optional (default `0`) + Split particles of the species when crossing the boundary from a lower + resolution domain to a higher resolution domain. + +* ``<species_name>.split_type`` (`int`) optional (default `0`) + Splitting technique. When `0`, particles are split along the simulation + axes (4 particles in 2D, 6 particles in 3D). When `1`, particles are split + along the diagonals (4 particles in 2D, 8 particles in 3D). + * ``warpx.serialize_ics`` (`0 or 1`) Whether or not to use OpenMP threading for particle initialization. @@ -392,6 +401,7 @@ Numerics and algorithms - ``0``: Vectorized version - ``1``: Non-optimized version + .. warning:: The vectorized version does not run on GPU. Use diff --git a/Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py b/Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py index 1006cdd77..169c56e7b 100755 --- a/Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py +++ b/Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py @@ -82,4 +82,4 @@ plt.tight_layout() plt.savefig('langmuir_multi_2d_analysis.png') # Automatically check the validity -assert overall_max_error < 0.02 +assert overall_max_error < 0.04 diff --git a/LICENSE.txt b/LICENSE.txt index 8135aac1b..1d49a0370 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,4 +1,4 @@ -WarpX v18.11 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.02 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 66aeff164..9116e49ec 100644 --- a/Python/setup.py +++ b/Python/setup.py @@ -20,7 +20,7 @@ else: package_data = {} setup (name = 'pywarpx', - version = '18.11', + version = '19.02', packages = ['pywarpx'], package_dir = {'pywarpx':'pywarpx'}, description = """Wrapper of WarpX""", diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 60647e89b..28b9ec225 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -192,6 +192,24 @@ particleTypes = electrons positrons analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_analysis.py analysisOutputImage = langmuir_multi_analysis.png +[Langmuir_multi_nodal] +buildDir = . +inputFile = Examples/Tests/Langmuir/inputs.multi.rt +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 4 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 1 +runtime_params = warpx.do_dynamic_scheduling=0 warpx.do_nodal=1 +particleTypes = electrons positrons +analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_analysis.py +analysisOutputImage = langmuir_multi_analysis.png + [Langmuir_multi_psatd] buildDir = . inputFile = Examples/Tests/Langmuir/inputs.multi.rt @@ -210,6 +228,24 @@ particleTypes = electrons positrons analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_analysis.py analysisOutputImage = langmuir_multi_analysis.png +[Langmuir_multi_2d_nodal] +buildDir = . +inputFile = Examples/Tests/Langmuir/inputs.multi.2d.rt +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 4 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +runtime_params = warpx.do_nodal=1 +particleTypes = electrons positrons +analysisRoutine = Examples/Tests/Langmuir/langmuir_multi_2d_analysis.py +analysisOutputImage = langmuir_multi_2d_analysis.png + [Langmuir_multi_2d_psatd] buildDir = . inputFile = Examples/Tests/Langmuir/inputs.multi.2d.rt @@ -354,59 +390,3 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons - -[UnitTest_CurrentDeposition] -buildDir = tests/CurrentDeposition -inputFile = inputs -dim = 3 -addToCompileString = -restartTest = 0 -useMPI = 1 -numprocs = 1 -useOMP = 0 -numthreads = 0 -compileTest = 0 -doVis = 0 -outputFile = plotfiles/plt00000 - -[UnitTest_FieldGather] -buildDir = tests/FieldGather -inputFile = inputs -dim = 3 -addToCompileString = -restartTest = 0 -useMPI = 1 -numprocs = 1 -useOMP = 0 -numthreads = 0 -compileTest = 0 -doVis = 0 -outputFile = plotfiles/plt00000 - -[UnitTest_FieldSolver] -buildDir = tests/FieldSolver -inputFile = inputs -dim = 3 -addToCompileString = -restartTest = 0 -useMPI = 1 -numprocs = 1 -useOMP = 0 -numthreads = 0 -compileTest = 0 -doVis = 0 -outputFile = plotfiles/plt00000 - -[UnitTest_ParticlePusher] -buildDir = tests/ParticlePusher -inputFile = inputs -dim = 3 -addToCompileString = -restartTest = 0 -useMPI = 1 -numprocs = 1 -useOMP = 0 -numthreads = 0 -compileTest = 0 -doVis = 0 -outputFile = plotfiles/plt00000 diff --git a/Regression/prepare_file_travis.py b/Regression/prepare_file_travis.py index 82b3c8e0e..725219677 100644 --- a/Regression/prepare_file_travis.py +++ b/Regression/prepare_file_travis.py @@ -18,6 +18,8 @@ with open('WarpX-tests.ini') as f: # Replace default folder name text = re.sub('/home/regtester/AMReX_RegTesting', test_dir, text) +# Remove the web directory +text = re.sub('[\w\-\/]*/web', '', text) # Add doComparison = 0 for each test text = re.sub( '\[(?P<name>.*)\]\nbuildDir = ', @@ -42,8 +44,6 @@ text = re.sub( '\[Python_Langmuir\]\n(.+\n)*', '', text) # Remove Langmuir_x/y/z test (too long; not that useful) text = re.sub( '\[Langmuir_[xyz]\]\n(.+\n)*', '', text) -# Skip unit tests (too long; not that useful) -text = re.sub( '\[UnitTest_[a-zA-Z]+\]\n(.+\n)*', '', text) # Remove tests that do not have the right dimension if dim is not None: diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package new file mode 100644 index 000000000..f27528c7d --- /dev/null +++ b/Source/BoundaryConditions/Make.package @@ -0,0 +1,6 @@ +CEXE_sources += PML.cpp WarpXEvolvePML.cpp +CEXE_headers += PML.H +F90EXE_sources += PML_routines.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/BoundaryConditions +VPATH_LOCATIONS += $(WARPX_HOME)/Source/BoundaryConditions diff --git a/Source/WarpXPML.H b/Source/BoundaryConditions/PML.H index 7a86bf0da..7a86bf0da 100644 --- a/Source/WarpXPML.H +++ b/Source/BoundaryConditions/PML.H diff --git a/Source/WarpXPML.cpp b/Source/BoundaryConditions/PML.cpp index b07d4164d..03373995d 100644 --- a/Source/WarpXPML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -1,5 +1,5 @@ -#include <WarpXPML.H> +#include <PML.H> #include <WarpX.H> #include <WarpXConst.H> diff --git a/Source/WarpX_pml.F90 b/Source/BoundaryConditions/PML_routines.F90 index 380e52934..380e52934 100644 --- a/Source/WarpX_pml.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp new file mode 100644 index 000000000..b0688b2c1 --- /dev/null +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -0,0 +1,81 @@ +#include <cmath> +#include <limits> + +#include <WarpX.H> +#include <WarpXConst.H> +#include <WarpX_f.H> +#ifdef WARPX_USE_PY +#include <WarpX_py.H> +#endif + +#ifdef BL_USE_SENSEI_INSITU +#include <AMReX_AmrMeshInSituBridge.H> +#endif + +using namespace amrex; + +void +WarpX::DampPML () +{ + for (int lev = 0; lev <= finest_level; ++lev) { + DampPML(lev); + } +} + +void +WarpX::DampPML (int lev) +{ + DampPML(lev, PatchType::fine); + if (lev > 0) DampPML(lev, PatchType::coarse); +} + +void +WarpX::DampPML (int lev, PatchType patch_type) +{ + if (!do_pml) return; + + BL_PROFILE("WarpX::DampPML()"); + + if (pml[lev]->ok()) + { + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() + : pml[lev]->GetMultiSigmaBox_cp(); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*pml_E[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + WRPX_DAMP_PML(tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), + WRPX_PML_TO_FORTRAN(sigba[mfi])); + + if (pml_F) { + const Box& tnd = mfi.nodaltilebox(); + WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), + BL_TO_FORTRAN_3D((*pml_F)[mfi]), + WRPX_PML_TO_FORTRAN(sigba[mfi])); + } + } + } +} diff --git a/Source/WarpXBoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index 452062137..17f2b161c 100644 --- a/Source/WarpXBoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -7,7 +7,7 @@ #include <AMReX_PlotFileUtil.H> #include <AMReX_ParallelDescriptor.H> -#include "ParticleContainer.H" +#include "MultiParticleContainer.H" #include "WarpXConst.H" /// diff --git a/Source/WarpXBoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 96903f992..6c3e3c58a 100644 --- a/Source/WarpXBoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -1,4 +1,4 @@ -#include "WarpXBoostedFrameDiagnostic.H" +#include "BoostedFrameDiagnostic.H" #include <AMReX_MultiFabUtil.H> #include "WarpX_f.H" #include "WarpX.H" diff --git a/Source/WarpX_boosted_frame.F90 b/Source/Diagnostics/BoostedFrame_module.F90 index 3bce1817e..3bce1817e 100644 --- a/Source/WarpX_boosted_frame.F90 +++ b/Source/Diagnostics/BoostedFrame_module.F90 diff --git a/Source/Diagnostics/ElectrostaticIO.cpp b/Source/Diagnostics/ElectrostaticIO.cpp new file mode 100644 index 000000000..8ed205b1a --- /dev/null +++ b/Source/Diagnostics/ElectrostaticIO.cpp @@ -0,0 +1,123 @@ +#include <AMReX_MGT_Solver.H> +#include <AMReX_stencil_types.H> + +#include <WarpX.H> +#include <WarpX_f.H> + +namespace +{ + const std::string level_prefix {"Level_"}; +} + +using namespace amrex; + +void +WarpX:: +WritePlotFileES (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, + const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi, + const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E) +{ + BL_PROFILE("WarpX::WritePlotFileES()"); + + VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); + VisMF::SetHeaderVersion(plotfile_headerversion); + + const std::string& plotfilename = amrex::Concatenate(plot_file,istep[0]); + + amrex::Print() << " Writing plotfile " << plotfilename << "\n"; + + const int nlevels = finestLevel()+1; + + { + Vector<std::string> varnames; + Vector<std::unique_ptr<MultiFab> > mf(finest_level+1); + + for (int lev = 0; lev <= finest_level; ++lev) { + int ncomp = 5; + const int ngrow = 0; + mf[lev].reset(new MultiFab(grids[lev], dmap[lev], ncomp, ngrow)); + + int dcomp = 0; + amrex::average_node_to_cellcenter(*mf[lev], dcomp, *rho[lev], 0, 1); + if (lev == 0) { + varnames.push_back("rho"); + } + dcomp += 1; + + amrex::average_node_to_cellcenter(*mf[lev], dcomp , *E[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *E[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *E[lev][0], 0, 1); + + if (lev == 0) { + varnames.push_back("Ex"); + varnames.push_back("Ey"); + varnames.push_back("Ez"); + } + dcomp += 3; + + amrex::average_node_to_cellcenter(*mf[lev], dcomp, *phi[lev], 0, 1); + if (lev == 0) { + varnames.push_back("phi"); + } + dcomp += 1; + } + + Vector<std::string> rfs(1,"raw_fields"); // pre-build raw_fields/ + amrex::WriteMultiLevelPlotfile(plotfilename, finest_level+1, + amrex::GetVecOfConstPtrs(mf), + varnames, Geom(), t_new[0], istep, refRatio(), + "HyperCLaw-V1.1", + "Level_", + "Cell", + rfs); + } + + { + const std::string raw_plotfilename = plotfilename + "/raw_fields"; + const int nlevels = finestLevel()+1; + for (int lev = 0; lev < nlevels; ++lev) { + const DistributionMapping& dm = DistributionMap(lev); + + MultiFab Ex( E[lev][0]->boxArray(), dm, 1, 0); + MultiFab Ey( E[lev][1]->boxArray(), dm, 1, 0); + MultiFab Ez( E[lev][2]->boxArray(), dm, 1, 0); + MultiFab charge_density(rho[lev]->boxArray(), dm, 1, 0); + MultiFab potential(phi[lev]->boxArray(), dm, 1, 0); + + MultiFab::Copy(Ex, *E[lev][0], 0, 0, 1, 0); + MultiFab::Copy(Ey, *E[lev][1], 0, 0, 1, 0); + MultiFab::Copy(Ez, *E[lev][2], 0, 0, 1, 0); + MultiFab::Copy(charge_density, *rho[lev], 0, 0, 1, 0); + MultiFab::Copy(potential, *phi[lev], 0, 0, 1, 0); + + VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex")); + VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey")); + VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez")); + VisMF::Write(charge_density, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "rho")); + VisMF::Write(potential, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "phi")); + } + } + + Vector<std::string> particle_varnames; + particle_varnames.push_back("weight"); + + particle_varnames.push_back("momentum_x"); + particle_varnames.push_back("momentum_y"); + particle_varnames.push_back("momentum_z"); + + particle_varnames.push_back("Ex"); + particle_varnames.push_back("Ey"); + particle_varnames.push_back("Ez"); + + particle_varnames.push_back("Bx"); + particle_varnames.push_back("By"); + particle_varnames.push_back("Bz"); + + mypc->Checkpoint(plotfilename, true, particle_varnames); + + WriteJobInfo(plotfilename); + + WriteWarpXHeader(plotfilename); + + VisMF::SetHeaderVersion(current_version); +} diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package new file mode 100644 index 000000000..b569666ff --- /dev/null +++ b/Source/Diagnostics/Make.package @@ -0,0 +1,9 @@ +CEXE_sources += WarpXIO.cpp +CEXE_sources += BoostedFrameDiagnostic.cpp +CEXE_sources += ParticleIO.cpp +CEXE_headers += BoostedFrameDiagnostic.H +CEXE_headers += ElectrostaticIO.cpp +F90EXE_sources += BoostedFrame_module.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics diff --git a/Source/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index ced2d6b9a..be9809bac 100644 --- a/Source/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -1,5 +1,5 @@ -#include <ParticleContainer.H> +#include <MultiParticleContainer.H> #include <WarpX.H> using namespace amrex; diff --git a/Source/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index e68b4a4e8..e68b4a4e8 100644 --- a/Source/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp diff --git a/Source/Evolve/Make.package b/Source/Evolve/Make.package new file mode 100644 index 000000000..adc114ccc --- /dev/null +++ b/Source/Evolve/Make.package @@ -0,0 +1,7 @@ +CEXE_sources += WarpXEvolveEM.cpp +ifeq ($(DO_ELECTROSTATIC),TRUE) + CEXE_sources += WarpXEvolveES.cpp +endif + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Evolve +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Evolve
\ No newline at end of file diff --git a/Source/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 6d917e36d..4c87d8321 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -16,22 +16,6 @@ using namespace amrex; void -WarpX::Evolve (int numsteps) { - BL_PROFILE_REGION("WarpX::Evolve()"); - -#ifdef WARPX_DO_ELECTROSTATIC - if (do_electrostatic) { - EvolveES(numsteps); - } else { - EvolveEM(numsteps); - } -#else - EvolveEM(numsteps); -#endif // WARPX_DO_ELECTROSTATIC - -} - -void WarpX::EvolveEM (int numsteps) { BL_PROFILE("WarpX::EvolveEM()"); @@ -448,429 +432,6 @@ WarpX::OneStep_sub1 (Real curtime) } void -WarpX::EvolveB (Real dt) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - EvolveB(lev, dt); - } -} - -void -WarpX::EvolveB (int lev, Real dt) -{ - BL_PROFILE("WarpX::EvolveB()"); - EvolveB(lev, PatchType::fine, dt); - if (lev > 0) - { - EvolveB(lev, PatchType::coarse, dt); - } -} - -void -WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) -{ - const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; - if (patch_type == PatchType::fine) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); - } - - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - - // Loop through the grids, and over the tiles within each grid -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - Real wt = amrex::second(); - - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - // Call picsar routine for each tile - warpx_push_bvec( - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &WarpX::maxwell_fdtd_solver_id); - - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (patch_type == PatchType::coarse) cbx.refine(rr); - wt = (amrex::second() - wt) / cbx.d_numPts();\ - FArrayBox* costfab = cost->fabPtr(mfi); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( cbx, work_box, - { - costfab->plus(wt, work_box); - }); - } - } - - if (do_pml && pml[lev]->ok()) - { - const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_B[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - WRPX_PUSH_PML_BVEC( - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &WarpX::maxwell_fdtd_solver_id); - } - } -} - -void -WarpX::EvolveE (Real dt) -{ - for (int lev = 0; lev <= finest_level; ++lev) - { - EvolveE(lev, dt); - } -} - -void -WarpX::EvolveE (int lev, Real dt) -{ - BL_PROFILE("WarpX::EvolveE()"); - EvolveE(lev, PatchType::fine, dt); - if (lev > 0) - { - EvolveE(lev, PatchType::coarse, dt); - } -} - -void -WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) -{ - const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; - const Real c2dt = (PhysConst::c*PhysConst::c) * dt; - - int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; - const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; - if (patch_type == PatchType::fine) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); - jx = current_fp[lev][0].get(); - jy = current_fp[lev][1].get(); - jz = current_fp[lev][2].get(); - F = F_fp[lev].get(); - } - else if (patch_type == PatchType::coarse) - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); - jx = current_cp[lev][0].get(); - jy = current_cp[lev][1].get(); - jz = current_cp[lev][2].get(); - F = F_cp[lev].get(); - } - - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - - // Loop through the grids, and over the tiles within each grid -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - Real wt = amrex::second(); - - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - - // Call picsar routine for each tile - warpx_push_evec( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - BL_TO_FORTRAN_3D((*jx)[mfi]), - BL_TO_FORTRAN_3D((*jy)[mfi]), - BL_TO_FORTRAN_3D((*jz)[mfi]), - &mu_c2_dt, - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (F) - { - warpx_push_evec_f( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*F)[mfi]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2], - &WarpX::maxwell_fdtd_solver_id); - } - - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (patch_type == PatchType::coarse) cbx.refine(rr); - wt = (amrex::second() - wt) / cbx.d_numPts(); - FArrayBox* costfab = cost->fabPtr(mfi); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( cbx, work_box, - { - costfab->plus(wt, work_box); - }); - } - } - - if (do_pml && pml[lev]->ok()) - { - if (F) pml[lev]->ExchangeF(patch_type, F); - - const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_E[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - - WRPX_PUSH_PML_EVEC( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (pml_F) - { - WRPX_PUSH_PML_EVEC_F( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_F )[mfi]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2], - &WarpX::maxwell_fdtd_solver_id); - } - } - } -} - -void -WarpX::EvolveF (Real dt, DtType dt_type) -{ - if (!do_dive_cleaning) return; - - for (int lev = 0; lev <= finest_level; ++lev) - { - EvolveF(lev, dt, dt_type); - } -} - -void -WarpX::EvolveF (int lev, Real dt, DtType dt_type) -{ - if (!do_dive_cleaning) return; - - EvolveF(lev, PatchType::fine, dt, dt_type); - if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); -} - -void -WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) -{ - if (!do_dive_cleaning) return; - - BL_PROFILE("WarpX::EvolveF()"); - - static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; - - int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; - const auto& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *rho, *F; - if (patch_type == PatchType::fine) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - rho = rho_fp[lev].get(); - F = F_fp[lev].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - rho = rho_cp[lev].get(); - F = F_cp[lev].get(); - } - - const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; - - MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); - ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); - MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); - - if (do_pml && pml[lev]->ok()) - { - const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_F, TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& bx = mfi.tilebox(); - WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_ANYD((*pml_F )[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2]); - } - } -} - -void -WarpX::DampPML () -{ - for (int lev = 0; lev <= finest_level; ++lev) { - DampPML(lev); - } -} - -void -WarpX::DampPML (int lev) -{ - DampPML(lev, PatchType::fine); - if (lev > 0) DampPML(lev, PatchType::coarse); -} - -void -WarpX::DampPML (int lev, PatchType patch_type) -{ - if (!do_pml) return; - - BL_PROFILE("WarpX::DampPML()"); - - if (pml[lev]->ok()) - { - const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() - : pml[lev]->GetMultiSigmaBox_cp(); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(*pml_E[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - WRPX_DAMP_PML(tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - WRPX_PML_TO_FORTRAN(sigba[mfi])); - - if (pml_F) { - const Box& tnd = mfi.nodaltilebox(); - WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), - BL_TO_FORTRAN_3D((*pml_F)[mfi]), - WRPX_PML_TO_FORTRAN(sigba[mfi])); - } - } - } -} - -void WarpX::PushParticlesandDepose (Real cur_time) { // Evolve particles to p^{n+1/2} and x^{n+1} diff --git a/Source/WarpXElectrostatic.cpp b/Source/Evolve/WarpXEvolveES.cpp index c4e37ce1a..e1e4a3929 100644 --- a/Source/WarpXElectrostatic.cpp +++ b/Source/Evolve/WarpXEvolveES.cpp @@ -343,114 +343,3 @@ void WarpX::computeE(Vector<std::array<std::unique_ptr<MultiFab>, 3> >& E, #endif } } - -void -WarpX:: -WritePlotFileES (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, - const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi, - const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E) -{ - BL_PROFILE("WarpX::WritePlotFileES()"); - - VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); - VisMF::SetHeaderVersion(plotfile_headerversion); - - const std::string& plotfilename = amrex::Concatenate(plot_file,istep[0]); - - amrex::Print() << " Writing plotfile " << plotfilename << "\n"; - - const int nlevels = finestLevel()+1; - - { - Vector<std::string> varnames; - Vector<std::unique_ptr<MultiFab> > mf(finest_level+1); - - for (int lev = 0; lev <= finest_level; ++lev) { - int ncomp = 5; - const int ngrow = 0; - mf[lev].reset(new MultiFab(grids[lev], dmap[lev], ncomp, ngrow)); - - int dcomp = 0; - amrex::average_node_to_cellcenter(*mf[lev], dcomp, *rho[lev], 0, 1); - if (lev == 0) { - varnames.push_back("rho"); - } - dcomp += 1; - - amrex::average_node_to_cellcenter(*mf[lev], dcomp , *E[lev][0], 0, 1); - amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *E[lev][0], 0, 1); - amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *E[lev][0], 0, 1); - - if (lev == 0) { - varnames.push_back("Ex"); - varnames.push_back("Ey"); - varnames.push_back("Ez"); - } - dcomp += 3; - - amrex::average_node_to_cellcenter(*mf[lev], dcomp, *phi[lev], 0, 1); - if (lev == 0) { - varnames.push_back("phi"); - } - dcomp += 1; - } - - Vector<std::string> rfs(1,"raw_fields"); // pre-build raw_fields/ - amrex::WriteMultiLevelPlotfile(plotfilename, finest_level+1, - amrex::GetVecOfConstPtrs(mf), - varnames, Geom(), t_new[0], istep, refRatio(), - "HyperCLaw-V1.1", - "Level_", - "Cell", - rfs); - } - - { - const std::string raw_plotfilename = plotfilename + "/raw_fields"; - const int nlevels = finestLevel()+1; - for (int lev = 0; lev < nlevels; ++lev) { - const DistributionMapping& dm = DistributionMap(lev); - - MultiFab Ex( E[lev][0]->boxArray(), dm, 1, 0); - MultiFab Ey( E[lev][1]->boxArray(), dm, 1, 0); - MultiFab Ez( E[lev][2]->boxArray(), dm, 1, 0); - MultiFab charge_density(rho[lev]->boxArray(), dm, 1, 0); - MultiFab potential(phi[lev]->boxArray(), dm, 1, 0); - - MultiFab::Copy(Ex, *E[lev][0], 0, 0, 1, 0); - MultiFab::Copy(Ey, *E[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Ez, *E[lev][2], 0, 0, 1, 0); - MultiFab::Copy(charge_density, *rho[lev], 0, 0, 1, 0); - MultiFab::Copy(potential, *phi[lev], 0, 0, 1, 0); - - VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex")); - VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey")); - VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez")); - VisMF::Write(charge_density, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "rho")); - VisMF::Write(potential, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "phi")); - } - } - - Vector<std::string> particle_varnames; - particle_varnames.push_back("weight"); - - particle_varnames.push_back("momentum_x"); - particle_varnames.push_back("momentum_y"); - particle_varnames.push_back("momentum_z"); - - particle_varnames.push_back("Ex"); - particle_varnames.push_back("Ey"); - particle_varnames.push_back("Ez"); - - particle_varnames.push_back("Bx"); - particle_varnames.push_back("By"); - particle_varnames.push_back("Bz"); - - mypc->Checkpoint(plotfilename, true, particle_varnames); - - WriteJobInfo(plotfilename); - - WriteWarpXHeader(plotfilename); - - VisMF::SetHeaderVersion(current_version); -} diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package new file mode 100644 index 000000000..42eb6a6f0 --- /dev/null +++ b/Source/FieldSolver/Make.package @@ -0,0 +1,15 @@ +CEXE_headers += WarpX_K.H +CEXE_sources += WarpXPushFieldsEM.cpp +ifeq ($(USE_PSATD),TRUE) + CEXE_sources += WarpXFFT.cpp + F90EXE_sources += WarpX_fft.F90 +endif +ifeq ($(USE_OPENBC_POISSON),TRUE) + F90EXE_sources += openbc_poisson_solver.F90 +endif +ifeq ($DO_ELECTROSTATIC,TRUE) + F90EXE_sources += solve_E_nodal.F90 +endif + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver diff --git a/Source/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 1fd6ffa1d..f09290f29 100644 --- a/Source/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -324,7 +324,7 @@ WarpX::InitFFTDataPlan (int lev) { warpx_fft_dataplan_init(&nox_fft, &noy_fft, &noz_fft, (*dataptr_fp_fft[lev])[mfi].data, &FFTData::N, - dx_fp.data(), &dt[lev], &fftw_plan_measure ); + dx_fp.data(), &dt[lev], &fftw_plan_measure, &WarpX::do_nodal ); } } else if (Efield_fp_fft[lev][0]->local_size() == 0) @@ -334,7 +334,8 @@ WarpX::InitFFTDataPlan (int lev) nullfftdata.reset(new FFTData()); warpx_fft_dataplan_init(&nox_fft, &noy_fft, &noz_fft, nullfftdata->data, &FFTData::N, - dx_fp.data(), &dt[lev], &fftw_plan_measure ); + dx_fp.data(), &dt[lev], &fftw_plan_measure, + &WarpX::do_nodal ); } else { diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp new file mode 100644 index 000000000..6d588ae86 --- /dev/null +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -0,0 +1,428 @@ + +#include <cmath> +#include <limits> + +#include <WarpX.H> +#include <WarpXConst.H> +#include <WarpX_f.H> +#include <WarpX_K.H> +#ifdef WARPX_USE_PY +#include <WarpX_py.H> +#endif + +#ifdef BL_USE_SENSEI_INSITU +#include <AMReX_AmrMeshInSituBridge.H> +#endif + +using namespace amrex; + +void +WarpX::EvolveB (Real dt) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + EvolveB(lev, dt); + } +} + +void +WarpX::EvolveB (int lev, Real dt) +{ + BL_PROFILE("WarpX::EvolveB()"); + EvolveB(lev, PatchType::fine, dt); + if (lev > 0) + { + EvolveB(lev, PatchType::coarse, dt); + } +} + +void +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +{ + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array<Real,3>& dx = WarpX::CellSize(patch_level); + Real dtsdx = dt/dx[0], dtsdy = dt/dx[1], dtsdz = dt/dx[2]; + + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; + if (patch_type == PatchType::fine) + { + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + } + + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Real wt = amrex::second(); + + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + if (do_nodal) { + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bx_nodal(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); + }); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_by_nodal(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); + }); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); + }); + } else { + // Call picsar routine for each tile + warpx_push_bvec( + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*Bx)[mfi]), + BL_TO_FORTRAN_3D((*By)[mfi]), + BL_TO_FORTRAN_3D((*Bz)[mfi]), + &dtsdx, &dtsdy, &dtsdz, + &WarpX::maxwell_fdtd_solver_id); + } + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (amrex::second() - wt) / cbx.d_numPts(); + auto costfab = cost->array(mfi); + amrex::ParallelFor(cbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + costfab(i,j,k) += wt; + }); + } + } + + if (do_pml && pml[lev]->ok()) + { + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*pml_B[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + WRPX_PUSH_PML_BVEC( + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), + &dtsdx, &dtsdy, &dtsdz, + &WarpX::maxwell_fdtd_solver_id); + } + } +} + +void +WarpX::EvolveE (Real dt) +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + EvolveE(lev, dt); + } +} + +void +WarpX::EvolveE (int lev, Real dt) +{ + BL_PROFILE("WarpX::EvolveE()"); + EvolveE(lev, PatchType::fine, dt); + if (lev > 0) + { + EvolveE(lev, PatchType::coarse, dt); + } +} + +void +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +{ + const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; + const Real c2dt = (PhysConst::c*PhysConst::c) * dt; + + int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array<Real,3>& dx = WarpX::CellSize(patch_level); + Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; + + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; + if (patch_type == PatchType::fine) + { + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + jx = current_fp[lev][0].get(); + jy = current_fp[lev][1].get(); + jz = current_fp[lev][2].get(); + F = F_fp[lev].get(); + } + else if (patch_type == PatchType::coarse) + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + jx = current_cp[lev][0].get(); + jy = current_cp[lev][1].get(); + jz = current_cp[lev][2].get(); + F = F_cp[lev].get(); + } + + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Real wt = amrex::second(); + + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + + if (do_nodal) { + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& jxfab = jx->array(mfi); + auto const& jyfab = jy->array(mfi); + auto const& jzfab = jz->array(mfi); + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_nodal(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_nodal(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); + }); + } else { + // Call picsar routine for each tile + warpx_push_evec( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*Bx)[mfi]), + BL_TO_FORTRAN_3D((*By)[mfi]), + BL_TO_FORTRAN_3D((*Bz)[mfi]), + BL_TO_FORTRAN_3D((*jx)[mfi]), + BL_TO_FORTRAN_3D((*jy)[mfi]), + BL_TO_FORTRAN_3D((*jz)[mfi]), + &mu_c2_dt, + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + } + + if (F) + { + warpx_push_evec_f( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*F)[mfi]), + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, + &WarpX::maxwell_fdtd_solver_id); + } + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (amrex::second() - wt) / cbx.d_numPts(); + auto costfab = cost->array(mfi); + amrex::ParallelFor(cbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + costfab(i,j,k) += wt; + }); + } + } + + if (do_pml && pml[lev]->ok()) + { + if (F) pml[lev]->ExchangeF(patch_type, F); + + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*pml_E[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + + WRPX_PUSH_PML_EVEC( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + + if (pml_F) + { + WRPX_PUSH_PML_EVEC_F( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_F )[mfi]), + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, + &WarpX::maxwell_fdtd_solver_id); + } + } + } +} + +void +WarpX::EvolveF (Real dt, DtType dt_type) +{ + if (!do_dive_cleaning) return; + + for (int lev = 0; lev <= finest_level; ++lev) + { + EvolveF(lev, dt, dt_type); + } +} + +void +WarpX::EvolveF (int lev, Real dt, DtType dt_type) +{ + if (!do_dive_cleaning) return; + + EvolveF(lev, PatchType::fine, dt, dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); +} + +void +WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +{ + if (!do_dive_cleaning) return; + + BL_PROFILE("WarpX::EvolveF()"); + + static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; + + int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const auto& dx = WarpX::CellSize(patch_level); + const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + + MultiFab *Ex, *Ey, *Ez, *rho, *F; + if (patch_type == PatchType::fine) + { + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + rho = rho_fp[lev].get(); + F = F_fp[lev].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + rho = rho_cp[lev].get(); + F = F_cp[lev].get(); + } + + const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; + + MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); + ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); + MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); + MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + + if (do_pml && pml[lev]->ok()) + { + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*pml_F, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& bx = mfi.tilebox(); + WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN_ANYD((*pml_F )[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]), + &dtsdx[0], &dtsdx[1], &dtsdx[2]); + } + } +} + diff --git a/Source/FieldSolver/WarpX_K.H b/Source/FieldSolver/WarpX_K.H new file mode 100644 index 000000000..2772b764d --- /dev/null +++ b/Source/FieldSolver/WarpX_K.H @@ -0,0 +1,97 @@ +#ifndef WARPX_K_H_ +#define WARPX_K_H_ + +#include <AMReX_FArrayBox.H> + +using namespace amrex; + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_bx_nodal (int j, int k, int l, Array4<Real> const& Bx, + Array4<Real const> const& Ey, Array4<Real const> const& Ez, + Real dtsdy, Real dtsdz) +{ +#if (AMREX_SPACEDIM == 3) + Bx(j,k,l) = Bx(j,k,l) - 0.5*dtsdy * (Ez(j,k+1,l ) - Ez(j,k-1,l )) + + 0.5*dtsdz * (Ey(j,k ,l+1) - Ey(j,k ,l-1)); +#else + Bx(j,k,0) = Bx(j,k,0) + 0.5*dtsdz * (Ey(j,k+1,0) - Ey(j,k-1,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_by_nodal (int j, int k, int l, Array4<Real> const& By, + Array4<Real const> const& Ex, Array4<Real const> const& Ez, + Real dtsdx, Real dtsdz) +{ +#if (AMREX_SPACEDIM == 3) + By(j,k,l) = By(j,k,l) + 0.5*dtsdx * (Ez(j+1,k,l ) - Ez(j-1,k,l )) + - 0.5*dtsdz * (Ex(j ,k,l+1) - Ex(j ,k,l-1)); +#else + By(j,k,0) = By(j,k,0) + 0.5*dtsdx * (Ez(j+1,k ,0) - Ez(j-1,k ,0)) + - 0.5*dtsdz * (Ex(j ,k+1,0) - Ex(j ,k-1,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_bz_nodal (int j, int k, int l, Array4<Real> const& Bz, + Array4<Real const> const& Ex, Array4<Real const> const& Ey, + Real dtsdx, Real dtsdy) +{ +#if (AMREX_SPACEDIM == 3) + Bz(j,k,l) = Bz(j,k,l) - 0.5*dtsdx * (Ey(j+1,k ,l) - Ey(j-1,k ,l)) + + 0.5*dtsdy * (Ex(j ,k+1,l) - Ex(j ,k-1,l)); +#else + Bz(j,k,0) = Bz(j,k,0) - 0.5*dtsdx * (Ey(j+1,k ,0) - Ey(j-1,k ,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ex_nodal (int j, int k, int l, Array4<Real> const& Ex, + Array4<Real const> const& By, Array4<Real const> const& Bz, + Array4<Real const> const& jx, + Real mudt, Real dtsdy, Real dtsdz) +{ +#if (AMREX_SPACEDIM == 3) + Ex(j,k,l) = Ex(j,k,l) + 0.5*dtsdy * (Bz(j,k+1,l ) - Bz(j,k-1,l )) + - 0.5*dtsdz * (By(j,k ,l+1) - By(j,k ,l-1)) + - mudt * jx(j,k,l); +#else + Ex(j,k,0) = Ex(j,k,0) - 0.5*dtsdz * (By(j,k+1,0) - By(j,k-1,0)) + - mudt * jx(j,k,0); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ey_nodal (int j, int k, int l, Array4<Real> const& Ey, + Array4<Real const> const& Bx, Array4<Real const> const& Bz, + Array4<Real const> const& jy, + Real mudt, Real dtsdx, Real dtsdz) +{ +#if (AMREX_SPACEDIM == 3) + Ey(j,k,l) = Ey(j,k,l) - 0.5*dtsdx * (Bz(j+1,k,l ) - Bz(j-1,k,l )) + + 0.5*dtsdz * (Bx(j ,k,l+1) - Bx(j ,k,l-1)) + - mudt * jy(j,k,l); +#else + Ey(j,k,0) = Ey(j,k,0) - 0.5*dtsdx * (Bz(j+1,k ,0) - Bz(j-1,k ,0)) + + 0.5*dtsdz * (Bx(j ,k+1,0) - Bx(j ,k-1,0)) + - mudt * jy(j,k,0); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ez_nodal (int j, int k, int l, Array4<Real> const& Ez, + Array4<Real const> const& Bx, Array4<Real const> const& By, + Array4<Real const> const& jz, + Real mudt, Real dtsdx, Real dtsdy) +{ +#if (AMREX_SPACEDIM == 3) + Ez(j,k,l) = Ez(j,k,l) + 0.5*dtsdx * (By(j+1,k ,l) - By(j-1,k ,l)) + - 0.5*dtsdy * (Bx(j ,k+1,l) - Bx(j ,k-1,l)) + - mudt * jz(j,k,l); +#else + Ez(j,k,0) = Ez(j,k,0) + 0.5*dtsdx * (By(j+1,k,0) - By(j-1,k,0)) + - mudt * jz(j,k,0); +#endif +} + +#endif diff --git a/Source/WarpX_fft.F90 b/Source/FieldSolver/WarpX_fft.F90 index e605b2be0..35357a63f 100644 --- a/Source/WarpX_fft.F90 +++ b/Source/FieldSolver/WarpX_fft.F90 @@ -81,7 +81,7 @@ contains !! !! Note: fft_data is a stuct containing 22 pointers to arrays !! 1-11: padded arrays in real space ; 12-22 arrays for the fields in Fourier space - subroutine warpx_fft_dataplan_init (nox, noy, noz, fft_data, ndata, dx_wrpx, dt_wrpx, fftw_measure) & + subroutine warpx_fft_dataplan_init (nox, noy, noz, fft_data, ndata, dx_wrpx, dt_wrpx, fftw_measure, do_nodal) & bind(c,name='warpx_fft_dataplan_init') USE picsar_precision, only: idp use shared_data, only : c_dim, p3dfft_flag, fftw_plan_measure, & @@ -104,6 +104,7 @@ contains integer, intent(in) :: nox, noy, noz, ndata integer, intent(in) :: fftw_measure + integer, intent(in) :: do_nodal type(c_ptr), intent(inout) :: fft_data(ndata) real(c_double), intent(in) :: dx_wrpx(3), dt_wrpx @@ -119,7 +120,11 @@ contains nzguards = 0_idp ! For the calculation of the modified [k] vectors - l_staggered = .TRUE. + if (do_nodal == 0) then + l_staggered = .TRUE. + else + l_staggered = .FALSE. + endif norderx = int(nox, idp) nordery = int(noy, idp) norderz = int(noz, idp) diff --git a/Source/openbc_f.F90 b/Source/FieldSolver/openbc_poisson_solver.F90 index c18b1db24..c18b1db24 100644 --- a/Source/openbc_f.F90 +++ b/Source/FieldSolver/openbc_poisson_solver.F90 diff --git a/Source/FieldSolver/solve_E_nodal.F90 b/Source/FieldSolver/solve_E_nodal.F90 new file mode 100644 index 000000000..e5ac50b2a --- /dev/null +++ b/Source/FieldSolver/solve_E_nodal.F90 @@ -0,0 +1,73 @@ +module warpx_ES_solve_E_nodal + + use iso_c_binding + use amrex_fort_module, only : amrex_real + + implicit none + +contains + +! This routine computes the node-centered electric field given a node-centered phi. +! The gradient is computed using 2nd-order centered differences. It assumes the +! Boundary conditions have already been set and that you have two rows of ghost cells +! for phi and one row of ghost cells for Ex, Ey, and Ez. +! Note that this routine includes the minus sign in E = - grad phi. +! +! Arguments: +! lo, hi: The corners of the valid box over which the gradient is taken +! Ex, Ey, Ez: The electric field in the x, y, and z directions. +! dx: The cell spacing +! + subroutine warpx_compute_E_nodal_3d (lo, hi, phi, Ex, Ey, Ez, dx) & + bind(c,name='warpx_compute_E_nodal_3d') + integer(c_int), intent(in) :: lo(3), hi(3) + real(amrex_real), intent(in) :: dx(3) + real(amrex_real), intent(in ) :: phi(lo(1)-2:hi(1)+2,lo(2)-2:hi(2)+2,lo(3)-2:hi(3)+2) + real(amrex_real), intent(inout) :: Ex (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) + real(amrex_real), intent(inout) :: Ey (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) + real(amrex_real), intent(inout) :: Ez (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) + + integer :: i, j, k + real(amrex_real) :: fac(3) + + fac = 0.5d0 / dx + + do k = lo(3)-1, hi(3)+1 + do j = lo(2)-1, hi(2)+1 + do i = lo(1)-1, hi(1)+1 + + Ex(i,j,k) = fac(1) * (phi(i-1,j,k) - phi(i+1,j,k)) + Ey(i,j,k) = fac(2) * (phi(i,j-1,k) - phi(i,j+1,k)) + Ez(i,j,k) = fac(3) * (phi(i,j,k-1) - phi(i,j,k+1)) + + end do + end do + end do + + end subroutine warpx_compute_E_nodal_3d + + subroutine warpx_compute_E_nodal_2d (lo, hi, phi, Ex, Ey, dx) & + bind(c,name='warpx_compute_E_nodal_2d') + integer(c_int), intent(in) :: lo(2), hi(2) + real(amrex_real), intent(in) :: dx(2) + real(amrex_real), intent(in ) :: phi(lo(1)-2:hi(1)+2,lo(2)-2:hi(2)+2) + real(amrex_real), intent(inout) :: Ex (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) + real(amrex_real), intent(inout) :: Ey (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) + + integer :: i, j + real(amrex_real) :: fac(2) + + fac = 0.5d0 / dx + + do j = lo(2)-1, hi(2)+1 + do i = lo(1)-1, hi(1)+1 + + Ex(i,j) = fac(1) * (phi(i-1,j) - phi(i+1,j)) + Ey(i,j) = fac(2) * (phi(i,j-1) - phi(i,j+1)) + + end do + end do + + end subroutine warpx_compute_E_nodal_2d + +end module warpx_ES_solve_E_nodal diff --git a/Source/Filter/Make.package b/Source/Filter/Make.package new file mode 100644 index 000000000..56c0e7852 --- /dev/null +++ b/Source/Filter/Make.package @@ -0,0 +1,4 @@ +F90EXE_sources += filter_module.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Filter +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Filter diff --git a/Source/WarpX_filter.F90 b/Source/Filter/filter_module.F90 index ea87be7ca..ea87be7ca 100644 --- a/Source/WarpX_filter.F90 +++ b/Source/Filter/filter_module.F90 diff --git a/Source/FortranInterface/Make.package b/Source/FortranInterface/Make.package new file mode 100644 index 000000000..fe0c2ba3d --- /dev/null +++ b/Source/FortranInterface/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += WarpX_f.H +F90EXE_sources += WarpX_f.F90 +F90EXE_sources += WarpX_picsar.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FortranInterface +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FortranInterface diff --git a/Source/WarpX_f.F90 b/Source/FortranInterface/WarpX_f.F90 index ccd6dd48a..ccd6dd48a 100644 --- a/Source/WarpX_f.F90 +++ b/Source/FortranInterface/WarpX_f.F90 diff --git a/Source/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index db7a5f9bd..1f5256483 100644 --- a/Source/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -476,7 +476,7 @@ extern "C" void warpx_fft_dataplan_init (const int* nox, const int* noy, const int* noz, void* fft_data, const int* ndata, const amrex_real* dx_w, const amrex_real* dt_w, - const int* fftw_plan_measure ); + const int* fftw_plan_measure, const int* do_nodal ); void warpx_fft_nullify (); void warpx_fft_push_eb (amrex_real* ex_w, const int* exlo, const int* exhi, amrex_real* ey_w, const int* eylo, const int* eyhi, diff --git a/Source/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 8e9d2db3a..8e9d2db3a 100644 --- a/Source/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 diff --git a/Source/CustomDensityProb.cpp b/Source/Initialization/CustomDensityProb.cpp index 3efcb13c5..3efcb13c5 100644 --- a/Source/CustomDensityProb.cpp +++ b/Source/Initialization/CustomDensityProb.cpp diff --git a/Source/CustomMomentumProb.cpp b/Source/Initialization/CustomMomentumProb.cpp index fa21252d0..fa21252d0 100644 --- a/Source/CustomMomentumProb.cpp +++ b/Source/Initialization/CustomMomentumProb.cpp diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package new file mode 100644 index 000000000..829b9a91c --- /dev/null +++ b/Source/Initialization/Make.package @@ -0,0 +1,9 @@ +CEXE_sources += CustomDensityProb.cpp +CEXE_sources += WarpXInitData.cpp +CEXE_sources += CustomMomentumProb.cpp +CEXE_sources += PlasmaInjector.cpp +CEXE_headers += PlasmaInjector.H +F90EXE_sources += WarpX_parser.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Initialization +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Initialization diff --git a/Source/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index c12607803..a25d1141b 100644 --- a/Source/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -239,6 +239,8 @@ public: using vec3 = std::array<amrex::Real, 3>; + PlasmaInjector(); + PlasmaInjector(int ispecies, const std::string& name); amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z); diff --git a/Source/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 832417686..3c120a7c1 100644 --- a/Source/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -201,6 +201,10 @@ void RegularPosition::getPositionUnitBox(vec3& r, int i_part, int ref_fac) r[2] = (0.5+iz_part)/nz; } +PlasmaInjector::PlasmaInjector(){ + part_pos = NULL; +} + PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) : species_id(ispecies), species_name(name) { diff --git a/Source/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 496b14e7a..ff5442b00 100644 --- a/Source/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -2,6 +2,7 @@ #include <numeric> #include <AMReX_ParallelDescriptor.H> +#include <AMReX_ParmParse.H> #include <WarpX.H> #include <WarpX_f.H> @@ -257,3 +258,77 @@ WarpX::InitOpenbc () } } #endif + +void +WarpX::InitLevelData (int lev, Real time) +{ + for (int i = 0; i < 3; ++i) { + current_fp[lev][i]->setVal(0.0); + Efield_fp[lev][i]->setVal(0.0); + Bfield_fp[lev][i]->setVal(0.0); + } + + if (lev > 0) { + for (int i = 0; i < 3; ++i) { + Efield_aux[lev][i]->setVal(0.0); + Bfield_aux[lev][i]->setVal(0.0); + + current_cp[lev][i]->setVal(0.0); + Efield_cp[lev][i]->setVal(0.0); + Bfield_cp[lev][i]->setVal(0.0); + } + } + + if (F_fp[lev]) { + F_fp[lev]->setVal(0.0); + } + + if (rho_fp[lev]) { + rho_fp[lev]->setVal(0.0); + } + + if (F_cp[lev]) { + F_cp[lev]->setVal(0.0); + } + + if (rho_cp[lev]) { + rho_cp[lev]->setVal(0.0); + } + + if (costs[lev]) { + costs[lev]->setVal(0.0); + } +} + +#ifdef WARPX_USE_PSATD + +void +WarpX::InitLevelDataFFT (int lev, Real time) +{ + Efield_fp_fft[lev][0]->setVal(0.0); + Efield_fp_fft[lev][1]->setVal(0.0); + Efield_fp_fft[lev][2]->setVal(0.0); + Bfield_fp_fft[lev][0]->setVal(0.0); + Bfield_fp_fft[lev][1]->setVal(0.0); + Bfield_fp_fft[lev][2]->setVal(0.0); + current_fp_fft[lev][0]->setVal(0.0); + current_fp_fft[lev][1]->setVal(0.0); + current_fp_fft[lev][2]->setVal(0.0); + rho_fp_fft[lev]->setVal(0.0); + + if (lev > 0) + { + Efield_cp_fft[lev][0]->setVal(0.0); + Efield_cp_fft[lev][1]->setVal(0.0); + Efield_cp_fft[lev][2]->setVal(0.0); + Bfield_cp_fft[lev][0]->setVal(0.0); + Bfield_cp_fft[lev][1]->setVal(0.0); + Bfield_cp_fft[lev][2]->setVal(0.0); + current_cp_fft[lev][0]->setVal(0.0); + current_cp_fft[lev][1]->setVal(0.0); + current_cp_fft[lev][2]->setVal(0.0); + rho_cp_fft[lev]->setVal(0.0); + } +} + +#endif diff --git a/Source/WarpX_parser.F90 b/Source/Initialization/WarpX_parser.F90 index 4b95c5ac7..4b95c5ac7 100644 --- a/Source/WarpX_parser.F90 +++ b/Source/Initialization/WarpX_parser.F90 diff --git a/Source/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 77a52dc26..77a52dc26 100644 --- a/Source/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H diff --git a/Source/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index cf669d32d..08d0ae861 100644 --- a/Source/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -7,7 +7,7 @@ #include <WarpX.H> #include <WarpXConst.H> #include <WarpX_f.H> -#include <ParticleContainer.H> +#include <MultiParticleContainer.H> using namespace amrex; @@ -35,10 +35,10 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies) std::transform(laser_type_s.begin(), laser_type_s.end(), laser_type_s.begin(), ::tolower); if (laser_type_s == "gaussian") { profile = laser_t::Gaussian; - } else if(laser_type_s == "harris") { - profile = laser_t::Harris; - } else if(laser_type_s == "parse_field_function") { - profile = laser_t::parse_field_function; + } else if(laser_type_s == "harris") { + profile = laser_t::Harris; + } else if(laser_type_s == "parse_field_function") { + profile = laser_t::parse_field_function; } else { amrex::Abort("Unknown laser type"); } diff --git a/Source/Laser/Make.package b/Source/Laser/Make.package new file mode 100644 index 000000000..96e1e8350 --- /dev/null +++ b/Source/Laser/Make.package @@ -0,0 +1,6 @@ +CEXE_sources += LaserParticleContainer.cpp +CEXE_headers += LaserParticleContainer.H +F90EXE_sources += WarpX_laser.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Laser +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Laser diff --git a/Source/WarpX_laser.F90 b/Source/Laser/WarpX_laser.F90 index b08413019..b08413019 100644 --- a/Source/WarpX_laser.F90 +++ b/Source/Laser/WarpX_laser.F90 diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 5d69dd0fa..75383e633 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -32,6 +32,18 @@ endif -include Make.package include $(WARPX_HOME)/Source/Make.package +include $(WARPX_HOME)/Source/BoundaryConditions/Make.package +include $(WARPX_HOME)/Source/Diagnostics/Make.package +include $(WARPX_HOME)/Source/FieldSolver/Make.package +include $(WARPX_HOME)/Source/Filter/Make.package +include $(WARPX_HOME)/Source/FortranInterface/Make.package +include $(WARPX_HOME)/Source/Initialization/Make.package +include $(WARPX_HOME)/Source/Laser/Make.package +include $(WARPX_HOME)/Source/Parallelization/Make.package +include $(WARPX_HOME)/Source/Particles/Make.package +include $(WARPX_HOME)/Source/Python/Make.package +include $(WARPX_HOME)/Source/Utils/Make.package +include $(WARPX_HOME)/Source/Evolve/Make.package include $(AMREX_HOME)/Src/Base/Make.package include $(AMREX_HOME)/Src/Particle/Make.package diff --git a/Source/Make.package b/Source/Make.package index f55e12297..3163a9de0 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -1,46 +1,9 @@ ifneq ($(USE_PYTHON_MAIN),TRUE) CEXE_sources += main.cpp -else - CEXE_sources += WarpXWrappers.cpp - CEXE_headers += WarpXWrappers.h endif -CEXE_sources += WarpX_py.cpp - -CEXE_sources += WarpX.cpp WarpXInitData.cpp WarpXEvolve.cpp WarpXIO.cpp WarpXProb.cpp WarpXRegrid.cpp -CEXE_sources += WarpXTagging.cpp WarpXComm.cpp WarpXMove.cpp WarpXBoostedFrameDiagnostic.cpp - -CEXE_sources += ParticleIO.cpp -CEXE_sources += ParticleContainer.cpp WarpXParticleContainer.cpp PhysicalParticleContainer.cpp LaserParticleContainer.cpp RigidInjectedParticleContainer.cpp - -CEXE_headers += WarpX_py.H - -CEXE_headers += WarpX.H WarpX_f.H WarpXConst.H WarpXBoostedFrameDiagnostic.H -CEXE_sources += WarpXConst.cpp - -CEXE_headers += ParticleContainer.H WarpXParticleContainer.H PhysicalParticleContainer.H LaserParticleContainer.H RigidInjectedParticleContainer.H - -CEXE_headers += PlasmaInjector.H -CEXE_sources += PlasmaInjector.cpp CustomDensityProb.cpp CustomMomentumProb.cpp - -CEXE_sources += WarpXPML.cpp WarpXUtil.cpp -CEXE_headers += WarpXPML.H WarpXUtil.H - -F90EXE_sources += WarpX_f.F90 WarpX_picsar.F90 WarpX_laser.F90 WarpX_pml.F90 WarpX_electrostatic.F90 -F90EXE_sources += WarpX_boosted_frame.F90 WarpX_filter.F90 WarpX_parser.F90 - -ifeq ($(USE_OPENBC_POISSON),TRUE) - F90EXE_sources += openbc_f.F90 -endif - -ifeq ($(USE_PSATD),TRUE) - CEXE_sources += WarpXFFT.cpp - F90EXE_sources += WarpX_fft.F90 -endif - -ifeq ($(DO_ELECTROSTATIC),TRUE) - CEXE_sources += WarpXElectrostatic.cpp -endif +CEXE_sources += WarpX.cpp +CEXE_headers += WarpX.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source VPATH_LOCATIONS += $(WARPX_HOME)/Source diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package new file mode 100644 index 000000000..cbb1b5234 --- /dev/null +++ b/Source/Parallelization/Make.package @@ -0,0 +1,5 @@ +CEXE_sources += WarpXComm.cpp +CEXE_sources += WarpXRegrid.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parallelization diff --git a/Source/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 69a57d2ac..09767c265 100644 --- a/Source/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -489,7 +489,7 @@ WarpX::SyncCurrent () const auto& period = Geom(lev).periodicity(); current_fp[lev][0]->OverrideSync(*current_fp_owner_masks[lev][0], period); current_fp[lev][1]->OverrideSync(*current_fp_owner_masks[lev][1], period); - current_fp[lev][2]->OverrideSync(*current_fp_owner_masks[lev][2],period); + current_fp[lev][2]->OverrideSync(*current_fp_owner_masks[lev][2], period); } for (int lev = 1; lev <= finest_level; ++lev) { diff --git a/Source/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 8d7873041..8d7873041 100644 --- a/Source/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package new file mode 100644 index 000000000..acbfe3390 --- /dev/null +++ b/Source/Particles/Make.package @@ -0,0 +1,14 @@ +F90EXE_sources += deposit_cic.F90 +F90EXE_sources += interpolate_cic.F90 +F90EXE_sources += push_particles_ES.F90 +CEXE_sources += PhysicalParticleContainer.cpp +CEXE_sources += MultiParticleContainer.cpp +CEXE_sources += WarpXParticleContainer.cpp +CEXE_sources += RigidInjectedParticleContainer.cpp +CEXE_headers += MultiParticleContainer.H +CEXE_headers += WarpXParticleContainer.H +CEXE_headers += RigidInjectedParticleContainer.H +CEXE_headers += PhysicalParticleContainer.H + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/ParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 9f732f355..b21247ff6 100644 --- a/Source/ParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -175,10 +175,13 @@ public: // Number of coefficients for the stencil of the NCI corrector. // The stencil is applied in the z direction only. static constexpr int nstencilz_fdtd_nci_corr=5; + amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_ex; amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_by; std::vector<std::string> GetSpeciesNames() const { return species_names; } + + PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } protected: @@ -195,6 +198,8 @@ private: // physical particles (+ laser) amrex::Vector<std::unique_ptr<WarpXParticleContainer> > allcontainers; + // Temporary particle container, used e.g. for particle splitting. + std::unique_ptr<PhysicalParticleContainer> pc_tmp; void ReadParameters (); diff --git a/Source/ParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 8aac692d6..1b644b543 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -2,7 +2,7 @@ #include <algorithm> #include <string> -#include <ParticleContainer.H> +#include <MultiParticleContainer.H> #include <WarpX_f.H> #include <WarpX.H> @@ -28,6 +28,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) if (WarpX::use_laser) { allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1)); } + pc_tmp.reset(new PhysicalParticleContainer(amr_core)); } void @@ -81,6 +82,7 @@ MultiParticleContainer::AllocData () for (auto& pc : allcontainers) { pc->AllocData(); } + pc_tmp->AllocData(); } void @@ -89,6 +91,7 @@ MultiParticleContainer::InitData () for (auto& pc : allcontainers) { pc->InitData(); } + pc_tmp->InitData(); } @@ -329,6 +332,7 @@ MultiParticleContainer::PostRestart () for (auto& pc : allcontainers) { pc->PostRestart(); } + pc_tmp->PostRestart(); } void diff --git a/Source/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index bd6015437..5847a6359 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -15,6 +15,9 @@ public: PhysicalParticleContainer (amrex::AmrCore* amr_core, int ispecies, const std::string& name); + + PhysicalParticleContainer (amrex::AmrCore* amr_core); + virtual ~PhysicalParticleContainer () {} virtual void InitData () override; @@ -77,6 +80,8 @@ public: virtual void PostRestart () final {} + void SplitParticles(int lev); + // Inject particles in Box 'part_box' virtual void AddParticles (int lev); void AddPlasma(int lev, amrex::RealBox part_realbox = amrex::RealBox()); diff --git a/Source/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 85644760c..b534de916 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1,7 +1,7 @@ #include <limits> #include <sstream> -#include <ParticleContainer.H> +#include <MultiParticleContainer.H> #include <WarpX_f.H> #include <WarpX.H> #include <WarpXConst.H> @@ -78,7 +78,15 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp ParmParse pp(species_name); pp.query("boost_adjust_transverse_positions", boost_adjust_transverse_positions); - pp.query("do_backward_propagation", do_backward_propagation); + pp.query("do_backward_propagation", do_backward_propagation); + pp.query("do_splitting", do_splitting); + pp.query("split_type", split_type); +} + +PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) + : WarpXParticleContainer(amr_core, 0) +{ + plasma_injector.reset(new PlasmaInjector()); } void PhysicalParticleContainer::InitData() @@ -1040,7 +1048,6 @@ PhysicalParticleContainer::Evolve (int lev, #else int thread_num = 0; #endif - if (local_rho[thread_num] == nullptr) local_rho[thread_num].reset( new amrex::FArrayBox()); if (local_jx[thread_num] == nullptr) local_jx[thread_num].reset( new amrex::FArrayBox()); if (local_jy[thread_num] == nullptr) local_jy[thread_num].reset( new amrex::FArrayBox()); @@ -1396,7 +1403,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_pxr_pp); PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], - m_giv[thread_num], dt); + m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_pxr_pp); // @@ -1426,6 +1433,162 @@ PhysicalParticleContainer::Evolve (int lev, } } } + // Split particles + if (do_splitting){ SplitParticles(lev); } +} + +// Loop over all particles in the particle container and +// split particles tagged with p.id()=DoSplitParticleID +void +PhysicalParticleContainer::SplitParticles(int lev) +{ + auto& mypc = WarpX::GetInstance().GetPartContainer(); + auto& pctmp_split = mypc.GetPCtmp(); + Cuda::DeviceVector<Real> xp, yp, zp; + RealVector psplit_x, psplit_y, psplit_z, psplit_w; + RealVector psplit_ux, psplit_uy, psplit_uz; + long np_split_to_add = 0; + long np_split; + if(split_type==0) + { + np_split = pow(2, AMREX_SPACEDIM); + } else { + np_split = 2*AMREX_SPACEDIM; + } + + // Loop over particle interator + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + pti.GetPosition(xp, yp, zp); + const std::array<Real,3>& dx = WarpX::CellSize(lev); + // particle Array Of Structs data + auto& particles = pti.GetArrayOfStructs(); + // particle Struct Of Arrays data + auto& attribs = pti.GetAttribs(); + auto& wp = attribs[PIdx::w ]; + auto& uxp = attribs[PIdx::ux]; + auto& uyp = attribs[PIdx::uy]; + auto& uzp = attribs[PIdx::uz]; + const long np = pti.numParticles(); + for(int i=0; i<np; i++){ + auto& p = particles[i]; + if (p.id() == DoSplitParticleID){ + // If particle is tagged, split it and put the + // split particles in local arrays psplit_x etc. + np_split_to_add += np_split; +#if (AMREX_SPACEDIM==2) + if (split_type==0){ + // Split particle in two along each axis + // 4 particles in 2d + for (int ishift = -1; ishift < 2; ishift +=2 ){ + for (int kshift = -1; kshift < 2; kshift +=2 ){ + // Add one particle with offset in x and z + psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_y.push_back( yp[i] ); + psplit_z.push_back( zp[i] + kshift*dx[2]/2 ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + } + } + } else { + // Split particle in two along each diagonal + // 4 particles in 2d + for (int ishift = -1; ishift < 2; ishift +=2 ){ + // Add one particle with offset in x + psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_y.push_back( yp[i] ); + psplit_z.push_back( zp[i] ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + // Add one particle with offset in z + psplit_x.push_back( xp[i] ); + psplit_y.push_back( yp[i] ); + psplit_z.push_back( zp[i] + ishift*dx[2]/2 ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + } + } +#elif (AMREX_SPACEDIM==3) + if (split_type==0){ + // Split particle in two along each axis + // 6 particles in 2d + for (int ishift = -1; ishift < 2; ishift +=2 ){ + for (int jshift = -1; jshift < 2; jshift +=2 ){ + for (int kshift = -1; kshift < 2; kshift +=2 ){ + // Add one particle with offset in x, y and z + psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_y.push_back( yp[i] + jshift*dx[1]/2 ); + psplit_z.push_back( zp[i] + jshift*dx[2]/2 ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + } + } + } + } else { + // Split particle in two along each diagonal + // 8 particles in 3d + for (int ishift = -1; ishift < 2; ishift +=2 ){ + // Add one particle with offset in x + psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_y.push_back( yp[i] ); + psplit_z.push_back( zp[i] ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + // Add one particle with offset in y + psplit_x.push_back( xp[i] ); + psplit_y.push_back( yp[i] + ishift*dx[1]/2 ); + psplit_z.push_back( zp[i] ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + // Add one particle with offset in z + psplit_x.push_back( xp[i] ); + psplit_y.push_back( yp[i] ); + psplit_z.push_back( zp[i] + ishift*dx[2]/2 ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + } + } +#endif + // invalidate the particle + p.m_idata.id = -p.m_idata.id; + } + } + } + // Add local arrays psplit_x etc. to the temporary + // particle container pctmp_split. Split particles + // are tagged with p.id()=NoSplitParticleID so that + // they are not re-split when entering a higher level + // AddNParticles calls Redistribute, so that particles + // in pctmp_split are in the proper grids and tiles + pctmp_split.AddNParticles(lev, + np_split_to_add, + psplit_x.dataPtr(), + psplit_y.dataPtr(), + psplit_z.dataPtr(), + psplit_ux.dataPtr(), + psplit_uy.dataPtr(), + psplit_uz.dataPtr(), + 1, + psplit_w.dataPtr(), + 1, NoSplitParticleID); + // Copy particles from tmp to current particle container + addParticles(pctmp_split,1); + // Clear tmp container + pctmp_split.clearParticles(); } void @@ -1543,6 +1706,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const int ll4symtry = false; long lvect_fieldgathe = 64; + warpx_geteb_energy_conserving( &np, m_xp[thread_num].dataPtr(), diff --git a/Source/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index 8a9ac9e6e..8a9ac9e6e 100644 --- a/Source/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H diff --git a/Source/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index db3623705..db3623705 100644 --- a/Source/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp diff --git a/Source/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 6b875c13f..d11ce49c9 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -184,7 +184,7 @@ public: void AddNParticles (int lev, int n, const amrex::Real* x, const amrex::Real* y, const amrex::Real* z, const amrex::Real* vx, const amrex::Real* vy, const amrex::Real* vz, - int nattr, const amrex::Real* attr, int uniqueparticles); + int nattr, const amrex::Real* attr, int uniqueparticles, int id=-1); void AddOneParticle (int lev, int grid, int tile, amrex::Real x, amrex::Real y, amrex::Real z, @@ -202,6 +202,11 @@ public: static int NextID () { return ParticleType::NextID(); } + bool do_splitting = false; + + // split along axes (0) or diagonals (1) + int split_type = 0; + protected: int species_id; @@ -219,7 +224,10 @@ protected: amrex::Vector<std::unique_ptr<amrex::FArrayBox> > local_jz; amrex::Vector<amrex::Cuda::DeviceVector<amrex::Real> > m_xp, m_yp, m_zp, m_giv; - + +private: + virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, + const int lev) override; }; #endif diff --git a/Source/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2bed6e0d3..1ba82e3d3 100644 --- a/Source/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1,7 +1,7 @@ #include <limits> -#include <ParticleContainer.H> +#include <MultiParticleContainer.H> #include <WarpXParticleContainer.H> #include <AMReX_AmrParGDB.H> #include <WarpX_f.H> @@ -98,7 +98,7 @@ WarpXParticleContainer::ReadParameters () #else do_tiling = true; #endif - pp.query("do_tiling", do_tiling); + pp.query("do_tiling", do_tiling); pp.query("do_not_push", do_not_push); initialized = true; @@ -152,8 +152,8 @@ void WarpXParticleContainer::AddNParticles (int lev, int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, - int nattr, const Real* attr, int uniqueparticles) -{ + int nattr, const Real* attr, int uniqueparticles, int id) +{ BL_ASSERT(nattr == 1); const Real* weight = attr; @@ -189,7 +189,12 @@ WarpXParticleContainer::AddNParticles (int lev, for (int i = ibegin; i < iend; ++i) { ParticleType p; - p.id() = ParticleType::NextID(); + if (id==-1) + { + p.id() = ParticleType::NextID(); + } else { + p.id() = id; + } p.cpu() = ParallelDescriptor::MyProc(); #if (AMREX_SPACEDIM == 3) p.pos(0) = x[i]; @@ -259,6 +264,8 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // WarpX assumes the same number of guard cells for Jx, Jy, Jz long ngJ = jx.nGrow(); + bool j_is_nodal = jx.is_nodal() and jy.is_nodal() and jz.is_nodal(); + // Deposit charge for particles that are not in the current buffers if (np_current > 0) { @@ -297,7 +304,68 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, auto jzntot = local_jz[thread_num]->length(); BL_PROFILE_VAR_START(blp_pxr_cd); - warpx_current_deposition( + if (j_is_nodal) { + const Real* p_wp = wp.dataPtr(); + const Real* p_gaminv = m_giv[thread_num].dataPtr(); + const Real* p_uxp = uxp.dataPtr(); + const Real* p_uyp = uyp.dataPtr(); + const Real* p_uzp = uzp.dataPtr(); + AsyncArray<Real> wptmp_aa(np_current); + Real* const wptmp = wptmp_aa.data(); + const Box& tile_box = pti.tilebox(); +#if (AMREX_SPACEDIM == 3) + const long nx = tile_box.length(0); + const long ny = tile_box.length(1); + const long nz = tile_box.length(2); +#else + const long nx = tile_box.length(0); + const long ny = 0; + const long nz = tile_box.length(1); +#endif + amrex::ParallelFor (np_current, [=] AMREX_GPU_DEVICE (long ip) { + wptmp[ip] = p_wp[ip] * p_gaminv[ip] * p_uxp[ip]; + }); + warpx_charge_deposition(jx_ptr, &np_current, + m_xp[thread_num].dataPtr(), + m_yp[thread_num].dataPtr(), + m_zp[thread_num].dataPtr(), + wptmp, + &this->charge, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, + &ngJ, &ngJ, &ngJ, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::current_deposition_algo); + amrex::ParallelFor (np_current, [=] AMREX_GPU_DEVICE (long ip) { + wptmp[ip] = p_wp[ip] * p_gaminv[ip] * p_uyp[ip]; + }); + warpx_charge_deposition(jy_ptr, &np_current, + m_xp[thread_num].dataPtr(), + m_yp[thread_num].dataPtr(), + m_zp[thread_num].dataPtr(), + wptmp, + &this->charge, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, + &ngJ, &ngJ, &ngJ, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::current_deposition_algo); + amrex::ParallelFor (np_current, [=] AMREX_GPU_DEVICE (long ip) { + wptmp[ip] = p_wp[ip] * p_gaminv[ip] * p_uzp[ip]; + }); + warpx_charge_deposition(jz_ptr, &np_current, + m_xp[thread_num].dataPtr(), + m_yp[thread_num].dataPtr(), + m_zp[thread_num].dataPtr(), + wptmp, + &this->charge, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, + &ngJ, &ngJ, &ngJ, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::current_deposition_algo); + } else { + warpx_current_deposition( jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), @@ -312,6 +380,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, &dt, &dx[0], &dx[1], &dx[2], &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect,&WarpX::current_deposition_algo); + } BL_PROFILE_VAR_STOP(blp_pxr_cd); BL_PROFILE_VAR_START(blp_accumulate); @@ -385,7 +454,68 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); - warpx_current_deposition( + if (j_is_nodal) { + const Real* p_wp = wp.dataPtr() + np_current; + const Real* p_gaminv = m_giv[thread_num].dataPtr() + np_current; + const Real* p_uxp = uxp.dataPtr() + np_current; + const Real* p_uyp = uyp.dataPtr() + np_current; + const Real* p_uzp = uzp.dataPtr() + np_current; + AsyncArray<Real> wptmp_aa(ncrse); + Real* const wptmp = wptmp_aa.data(); + const Box& tile_box = pti.tilebox(); +#if (AMREX_SPACEDIM == 3) + const long nx = tile_box.length(0); + const long ny = tile_box.length(1); + const long nz = tile_box.length(2); +#else + const long nx = tile_box.length(0); + const long ny = 0; + const long nz = tile_box.length(1); +#endif + amrex::ParallelFor (ncrse, [=] AMREX_GPU_DEVICE (long ip) { + wptmp[ip] = p_wp[ip] * p_gaminv[ip] * p_uxp[ip]; + }); + warpx_charge_deposition(jx_ptr, &ncrse, + m_xp[thread_num].dataPtr() +np_current, + m_yp[thread_num].dataPtr() +np_current, + m_zp[thread_num].dataPtr() +np_current, + wptmp, + &this->charge, + &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], + &cdx[0], &cdx[1], &cdx[2], &nx, &ny, &nz, + &ngJ, &ngJ, &ngJ, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::current_deposition_algo); + amrex::ParallelFor (ncrse, [=] AMREX_GPU_DEVICE (long ip) { + wptmp[ip] = p_wp[ip] * p_gaminv[ip] * p_uyp[ip]; + }); + warpx_charge_deposition(jy_ptr, &ncrse, + m_xp[thread_num].dataPtr() +np_current, + m_yp[thread_num].dataPtr() +np_current, + m_zp[thread_num].dataPtr() +np_current, + wptmp, + &this->charge, + &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], + &cdx[0], &cdx[1], &cdx[2], &nx, &ny, &nz, + &ngJ, &ngJ, &ngJ, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::current_deposition_algo); + amrex::ParallelFor (ncrse, [=] AMREX_GPU_DEVICE (long ip) { + wptmp[ip] = p_wp[ip] * p_gaminv[ip] * p_uzp[ip]; + }); + warpx_charge_deposition(jz_ptr, &ncrse, + m_xp[thread_num].dataPtr() +np_current, + m_yp[thread_num].dataPtr() +np_current, + m_zp[thread_num].dataPtr() +np_current, + wptmp, + &this->charge, + &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], + &cdx[0], &cdx[1], &cdx[2], &nx, &ny, &nz, + &ngJ, &ngJ, &ngJ, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::current_deposition_algo); + } else { + warpx_current_deposition( jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), @@ -402,6 +532,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, &dt, &cdx[0], &cdx[1], &cdx[2], &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect,&WarpX::current_deposition_algo); + } BL_PROFILE_VAR_STOP(blp_pxr_cd); BL_PROFILE_VAR_START(blp_accumulate); @@ -919,3 +1050,23 @@ WarpXParticleContainer::PushX (int lev, Real dt) } } } + +// This function is called in Redistribute, just after locate +void +WarpXParticleContainer::particlePostLocate(ParticleType& p, + const ParticleLocData& pld, + const int lev) +{ + // Tag particle if goes to higher level. + // It will be split later in the loop + if (pld.m_lev == lev+1 + and p.m_idata.id != NoSplitParticleID + and p.m_idata.id >= 0) + { + p.m_idata.id = DoSplitParticleID; + } + // For the moment, do not do anything if particles goes + // to lower level. + if (pld.m_lev == lev-1){ + } +} diff --git a/Source/Particles/deposit_cic.F90 b/Source/Particles/deposit_cic.F90 new file mode 100644 index 000000000..e4e1ace03 --- /dev/null +++ b/Source/Particles/deposit_cic.F90 @@ -0,0 +1,134 @@ +module warpx_ES_deposit_cic + + use iso_c_binding + use amrex_fort_module, only : amrex_real + + implicit none + +contains + +! This routine computes the charge density due to the particles using cloud-in-cell +! deposition. The Fab rho is assumed to be node-centered. +! +! Arguments: +! particles : a pointer to the particle array-of-structs +! ns : the stride length of particle struct (the size of the struct in number of reals) +! np : the number of particles +! weights : the particle weights +! charge : the charge of this particle species +! rho : a Fab that will contain the charge density on exit +! lo : a pointer to the lo corner of this valid box for rho, in index space +! hi : a pointer to the hi corner of this valid box for rho, in index space +! plo : the real position of the left-hand corner of the problem domain +! dx : the mesh spacing +! ng : the number of ghost cells in rho +! + subroutine warpx_deposit_cic_3d(particles, ns, np, & + weights, charge, rho, lo, hi, plo, dx, & + ng) & + bind(c,name='warpx_deposit_cic_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(in) :: weights(np) + real(amrex_real), intent(in) :: charge + integer, intent(in) :: lo(3) + integer, intent(in) :: hi(3) + integer, intent(in) :: ng + real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: plo(3) + real(amrex_real), intent(in) :: dx(3) + + integer i, j, k, n + real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi + real(amrex_real) lx, ly, lz + real(amrex_real) inv_dx(3) + real(amrex_real) qp, inv_vol + + inv_dx = 1.0d0/dx + + inv_vol = inv_dx(1) * inv_dx(2) * inv_dx(3) + + do n = 1, np + + qp = weights(n) * charge * inv_vol + + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + lz = (particles(3, n) - plo(3))*inv_dx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + rho(i, j, k ) = rho(i, j, k ) + wx_lo*wy_lo*wz_lo*qp + rho(i, j, k+1) = rho(i, j, k+1) + wx_lo*wy_lo*wz_hi*qp + rho(i, j+1, k ) = rho(i, j+1, k ) + wx_lo*wy_hi*wz_lo*qp + rho(i, j+1, k+1) = rho(i, j+1, k+1) + wx_lo*wy_hi*wz_hi*qp + rho(i+1, j, k ) = rho(i+1, j, k ) + wx_hi*wy_lo*wz_lo*qp + rho(i+1, j, k+1) = rho(i+1, j, k+1) + wx_hi*wy_lo*wz_hi*qp + rho(i+1, j+1, k ) = rho(i+1, j+1, k ) + wx_hi*wy_hi*wz_lo*qp + rho(i+1, j+1, k+1) = rho(i+1, j+1, k+1) + wx_hi*wy_hi*wz_hi*qp + + end do + + end subroutine warpx_deposit_cic_3d + + subroutine warpx_deposit_cic_2d(particles, ns, np, & + weights, charge, rho, lo, hi, plo, dx, & + ng) & + bind(c,name='warpx_deposit_cic_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(in) :: weights(np) + real(amrex_real), intent(in) :: charge + integer, intent(in) :: lo(2) + integer, intent(in) :: hi(2) + integer, intent(in) :: ng + real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: plo(2) + real(amrex_real), intent(in) :: dx(2) + + integer i, j, n + real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi + real(amrex_real) lx, ly + real(amrex_real) inv_dx(2) + real(amrex_real) qp, inv_vol + + inv_dx = 1.0d0/dx + + inv_vol = inv_dx(1) * inv_dx(2) + + do n = 1, np + + qp = weights(n) * charge * inv_vol + + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + + i = floor(lx) + j = floor(ly) + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + rho(i, j ) = rho(i, j ) + wx_lo*wy_lo*qp + rho(i, j+1) = rho(i, j+1) + wx_lo*wy_hi*qp + rho(i+1, j ) = rho(i+1, j ) + wx_hi*wy_lo*qp + rho(i+1, j+1) = rho(i+1, j+1) + wx_hi*wy_hi*qp + + end do + + end subroutine warpx_deposit_cic_2d + +end module warpx_ES_deposit_cic diff --git a/Source/Particles/interpolate_cic.F90 b/Source/Particles/interpolate_cic.F90 new file mode 100644 index 000000000..bc2c4f37e --- /dev/null +++ b/Source/Particles/interpolate_cic.F90 @@ -0,0 +1,371 @@ +module warpx_ES_interpolate_cic + + use iso_c_binding + use amrex_fort_module, only : amrex_real + + implicit none + +contains + +! This routine interpolates the electric field to the particle positions +! using cloud-in-cell interpolation. The electric fields are assumed to be +! node-centered. +! +! Arguments: +! particles : a pointer to the particle array-of-structs +! ns : the stride length of particle struct (the size of the struct in number of reals) +! np : the number of particles +! Ex_p : the electric field in the x-direction at the particle positions (output) +! Ey_p : the electric field in the y-direction at the particle positions (output) +! Ez_p : the electric field in the z-direction at the particle positions (output) +! Ex, Ey, Ez: Fabs conting the electric field on the mesh +! lo : a pointer to the lo corner of this valid box, in index space +! hi : a pointer to the hi corner of this valid box, in index space +! plo : the real position of the left-hand corner of the problem domain +! dx : the mesh spacing +! ng : the number of ghost cells for the E-field +! + subroutine warpx_interpolate_cic_3d(particles, ns, np, & + Ex_p, Ey_p, Ez_p, & + Ex, Ey, Ez, & + lo, hi, plo, dx, ng) & + bind(c,name='warpx_interpolate_cic_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) + integer, intent(in) :: ng + integer, intent(in) :: lo(3) + integer, intent(in) :: hi(3) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ez(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: plo(3) + real(amrex_real), intent(in) :: dx(3) + + integer i, j, k, n + real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi + real(amrex_real) lx, ly, lz + real(amrex_real) inv_dx(3) + inv_dx = 1.0d0/dx + + do n = 1, np + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + lz = (particles(3, n) - plo(3))*inv_dx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + Ex_p(n) = wx_lo*wy_lo*wz_lo*Ex(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ex(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ex(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ex(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ex(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) + + Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ey(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ey(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) + + Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ez(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ez(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) + + end do + + end subroutine warpx_interpolate_cic_3d + + + subroutine warpx_interpolate_cic_2d(particles, ns, np, & + Ex_p, Ey_p, & + Ex, Ey, & + lo, hi, plo, dx, ng) & + bind(c,name='warpx_interpolate_cic_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) + integer, intent(in) :: ng + integer, intent(in) :: lo(2) + integer, intent(in) :: hi(2) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: plo(2) + real(amrex_real), intent(in) :: dx(2) + + integer i, j, n + real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi + real(amrex_real) lx, ly + real(amrex_real) inv_dx(2) + inv_dx = 1.0d0/dx + + do n = 1, np + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + + i = floor(lx) + j = floor(ly) + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + Ex_p(n) = wx_lo*wy_lo*Ex(i, j ) + & + wx_lo*wy_hi*Ex(i, j+1) + & + wx_hi*wy_lo*Ex(i+1, j ) + & + wx_hi*wy_hi*Ex(i+1, j+1) + + Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & + wx_lo*wy_hi*Ey(i, j+1) + & + wx_hi*wy_lo*Ey(i+1, j ) + & + wx_hi*wy_hi*Ey(i+1, j+1) + + end do + + end subroutine warpx_interpolate_cic_2d + + + subroutine warpx_interpolate_cic_two_levels_3d(particles, ns, np, & + Ex_p, Ey_p, Ez_p, & + Ex, Ey, Ez, & + lo, hi, dx, & + cEx, cEy, cEz, & + mask, & + clo, chi, cdx, & + plo, ng, lev) & + bind(c,name='warpx_interpolate_cic_two_levels_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) + integer, intent(in) :: ng, lev + integer, intent(in) :: lo(3), hi(3) + integer, intent(in) :: clo(3), chi(3) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: Ez(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) + real(amrex_real), intent(in) :: cEx(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) + real(amrex_real), intent(in) :: cEy(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) + real(amrex_real), intent(in) :: cEz(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) + integer(c_int), intent(in) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) + real(amrex_real), intent(in) :: plo(3) + real(amrex_real), intent(in) :: dx(3), cdx(3) + + integer i, j, k, n + real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi + real(amrex_real) lx, ly, lz + real(amrex_real) inv_dx(3), inv_cdx(3) + inv_dx = 1.0d0/dx + inv_cdx = 1.0d0/cdx + + do n = 1, np + + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + lz = (particles(3, n) - plo(3))*inv_dx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + +! use the coarse E if near the level boundary + if (lev .eq. 1 .and. mask(i,j,k) .eq. 1) then + + lx = (particles(1, n) - plo(1))*inv_cdx(1) + ly = (particles(2, n) - plo(2))*inv_cdx(2) + lz = (particles(3, n) - plo(3))*inv_cdx(3) + + i = floor(lx) + j = floor(ly) + k = floor(lz) + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + Ex_p(n) = wx_lo*wy_lo*wz_lo*cEx(i, j, k ) + & + wx_lo*wy_lo*wz_hi*cEx(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*cEx(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*cEx(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*cEx(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*cEx(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*cEx(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*cEx(i+1, j+1, k+1) + + Ey_p(n) = wx_lo*wy_lo*wz_lo*cEy(i, j, k ) + & + wx_lo*wy_lo*wz_hi*cEy(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*cEy(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*cEy(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*cEy(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*cEy(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*cEy(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*cEy(i+1, j+1, k+1) + + Ez_p(n) = wx_lo*wy_lo*wz_lo*cEz(i, j, k ) + & + wx_lo*wy_lo*wz_hi*cEz(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*cEz(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*cEz(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*cEz(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*cEz(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*cEz(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*cEz(i+1, j+1, k+1) + +! otherwise use the fine + else + + wx_hi = lx - i + wy_hi = ly - j + wz_hi = lz - k + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + wz_lo = 1.0d0 - wz_hi + + Ex_p(n) = wx_lo*wy_lo*wz_lo*Ex(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ex(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ex(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ex(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ex(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) + + Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ey(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ey(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) + + Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & + wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & + wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & + wx_lo*wy_hi*wz_hi*Ez(i, j+1, k+1) + & + wx_hi*wy_lo*wz_lo*Ez(i+1, j, k ) + & + wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & + wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & + wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) + + end if + + end do + + end subroutine warpx_interpolate_cic_two_levels_3d + + + subroutine warpx_interpolate_cic_two_levels_2d(particles, ns, np, & + Ex_p, Ey_p, & + Ex, Ey, & + lo, hi, dx, & + cEx, cEy, & + mask, & + clo, chi, cdx, & + plo, ng, lev) & + bind(c,name='warpx_interpolate_cic_two_levels_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) + integer, intent(in) :: ng, lev + integer, intent(in) :: lo(2), hi(2) + integer, intent(in) :: clo(2), chi(2) + real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) + real(amrex_real), intent(in) :: cEx(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng) + real(amrex_real), intent(in) :: cEy(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng) + integer(c_int), intent(in) :: mask (lo(1):hi(1),lo(2):hi(2)) + real(amrex_real), intent(in) :: plo(2) + real(amrex_real), intent(in) :: dx(2), cdx(2) + + integer i, j, n + real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi + real(amrex_real) lx, ly + real(amrex_real) inv_dx(2), inv_cdx(2) + inv_dx = 1.0d0/dx + inv_cdx = 1.0d0/cdx + + do n = 1, np + + lx = (particles(1, n) - plo(1))*inv_dx(1) + ly = (particles(2, n) - plo(2))*inv_dx(2) + + i = floor(lx) + j = floor(ly) + +! use the coarse E if near the level boundary + if (lev .eq. 1 .and. mask(i,j) .eq. 1) then + + lx = (particles(1, n) - plo(1))*inv_cdx(1) + ly = (particles(2, n) - plo(2))*inv_cdx(2) + + i = floor(lx) + j = floor(ly) + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + Ex_p(n) = wx_lo*wy_lo*cEx(i, j ) + & + wx_lo*wy_hi*cEx(i, j+1) + & + wx_hi*wy_lo*cEx(i+1, j ) + & + wx_hi*wy_hi*cEx(i+1, j+1) + + Ey_p(n) = wx_lo*wy_lo*cEy(i, j ) + & + wx_lo*wy_hi*cEy(i, j+1) + & + wx_hi*wy_lo*cEy(i+1, j ) + & + wx_hi*wy_hi*cEy(i+1, j+1) + +! otherwise use the fine + else + + wx_hi = lx - i + wy_hi = ly - j + + wx_lo = 1.0d0 - wx_hi + wy_lo = 1.0d0 - wy_hi + + Ex_p(n) = wx_lo*wy_lo*Ex(i, j ) + & + wx_lo*wy_hi*Ex(i, j+1) + & + wx_hi*wy_lo*Ex(i+1, j ) + & + wx_hi*wy_hi*Ex(i+1, j+1) + + Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & + wx_lo*wy_hi*Ey(i, j+1) + & + wx_hi*wy_lo*Ey(i+1, j ) + & + wx_hi*wy_hi*Ey(i+1, j+1) + + end if + + end do + + end subroutine warpx_interpolate_cic_two_levels_2d + +end module warpx_ES_interpolate_cic diff --git a/Source/Particles/push_particles_ES.F90 b/Source/Particles/push_particles_ES.F90 new file mode 100644 index 000000000..9c7bd9e92 --- /dev/null +++ b/Source/Particles/push_particles_ES.F90 @@ -0,0 +1,258 @@ +module warpx_ES_push_particles + + use iso_c_binding + use amrex_fort_module, only : amrex_real + + implicit none + +contains + +! +! This routine updates the particle positions and velocities using the +! leapfrog time integration algorithm, given the electric fields at the +! particle positions. It also enforces specular reflection off the domain +! walls. +! +! Arguments: +! particles : a pointer to the particle array-of-structs +! ns : the stride length of particle struct (the size of the struct in number of reals) +! np : the number of particles +! vx_p : the particle x-velocities +! vy_p : the particle y-velocities +! vz_p : the particle z-velocities +! Ex_p : the electric field in the x-direction at the particle positions +! Ey_p : the electric field in the y-direction at the particle positions +! Ez_p : the electric field in the z-direction at the particle positions +! charge : the charge of this particle species +! mass : the mass of this particle species +! dt : the time step +! prob_lo : the left-hand corner of the problem domain +! prob_hi : the right-hand corner of the problem domain +! + subroutine warpx_push_leapfrog_3d(particles, ns, np, & + vx_p, vy_p, vz_p, & + Ex_p, Ey_p, Ez_p, & + charge, mass, dt, & + prob_lo, prob_hi) & + bind(c,name='warpx_push_leapfrog_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(inout) :: particles(ns,np) + real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) + real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np), Ez_p(np) + real(amrex_real), intent(in) :: charge + real(amrex_real), intent(in) :: mass + real(amrex_real), intent(in) :: dt + real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) + + integer n + real(amrex_real) fac + + fac = charge * dt / mass + + do n = 1, np + + vx_p(n) = vx_p(n) + fac * Ex_p(n) + vy_p(n) = vy_p(n) + fac * Ey_p(n) + vz_p(n) = vz_p(n) + fac * Ez_p(n) + + particles(1, n) = particles(1, n) + dt * vx_p(n) + particles(2, n) = particles(2, n) + dt * vy_p(n) + particles(3, n) = particles(3, n) + dt * vz_p(n) + +! bounce off the walls in the x... + do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) + if (particles(1, n) .lt. prob_lo(1)) then + particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) + else + particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) + end if + vx_p(n) = -vx_p(n) + end do + +! ... y... + do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) + if (particles(2, n) .lt. prob_lo(2)) then + particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) + else + particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) + end if + vy_p(n) = -vy_p(n) + end do + +! ... and z directions + do while (particles(3, n) .lt. prob_lo(3) .or. particles(3, n) .gt. prob_hi(3)) + if (particles(3, n) .lt. prob_lo(3)) then + particles(3, n) = 2.d0*prob_lo(3) - particles(3, n) + else + particles(3, n) = 2.d0*prob_hi(3) - particles(3, n) + end if + vz_p(n) = -vz_p(n) + end do + + end do + + end subroutine warpx_push_leapfrog_3d + + subroutine warpx_push_leapfrog_2d(particles, ns, np, & + vx_p, vy_p, & + Ex_p, Ey_p, & + charge, mass, dt, & + prob_lo, prob_hi) & + bind(c,name='warpx_push_leapfrog_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(inout) :: particles(ns,np) + real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) + real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np) + real(amrex_real), intent(in) :: charge + real(amrex_real), intent(in) :: mass + real(amrex_real), intent(in) :: dt + real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) + + integer n + real(amrex_real) fac + + fac = charge * dt / mass + + do n = 1, np + + vx_p(n) = vx_p(n) + fac * Ex_p(n) + vy_p(n) = vy_p(n) + fac * Ey_p(n) + + particles(1, n) = particles(1, n) + dt * vx_p(n) + particles(2, n) = particles(2, n) + dt * vy_p(n) + +! bounce off the walls in the x... + do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) + if (particles(1, n) .lt. prob_lo(1)) then + particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) + else + particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) + end if + vx_p(n) = -vx_p(n) + end do + +! ... y... + do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) + if (particles(2, n) .lt. prob_lo(2)) then + particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) + else + particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) + end if + vy_p(n) = -vy_p(n) + end do + + end do + + end subroutine warpx_push_leapfrog_2d + + +! +! This routine advances the particle positions using the current +! velocity. This is needed to desynchronize the particle positions +! from the velocities after particle initialization. +! +! Arguments: +! particles : a pointer to the particle array-of-structs +! ns : the stride length of particle struct (the size of the struct in number of reals) +! np : the number of particles +! xx_p : the electric field in the x-direction at the particle positions +! vy_p : the electric field in the y-direction at the particle positions +! vz_p : the electric field in the z-direction at the particle positions +! dt : the time step +! prob_lo : the left-hand corner of the problem domain +! prob_hi : the right-hand corner of the problem domain +! + subroutine warpx_push_leapfrog_positions_3d(particles, ns, np, & + vx_p, vy_p, vz_p, dt, & + prob_lo, prob_hi) & + bind(c,name='warpx_push_leapfrog_positions_3d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(inout) :: particles(ns,np) + real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) + real(amrex_real), intent(in) :: dt + real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) + + integer n + + do n = 1, np + + particles(1, n) = particles(1, n) + dt * vx_p(n) + particles(2, n) = particles(2, n) + dt * vy_p(n) + particles(3, n) = particles(3, n) + dt * vz_p(n) + +! bounce off the walls in the x... + do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) + if (particles(1, n) .lt. prob_lo(1)) then + particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) + else + particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) + end if + vx_p(n) = -vx_p(n) + end do + +! ... y... + do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) + if (particles(2, n) .lt. prob_lo(2)) then + particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) + else + particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) + end if + vy_p(n) = -vy_p(n) + end do + +! ... and z directions + do while (particles(3, n) .lt. prob_lo(3) .or. particles(3, n) .gt. prob_hi(3)) + if (particles(3, n) .lt. prob_lo(3)) then + particles(3, n) = 2.d0*prob_lo(3) - particles(3, n) + else + particles(3, n) = 2.d0*prob_hi(3) - particles(3, n) + end if + vz_p(n) = -vz_p(n) + end do + + end do + + end subroutine warpx_push_leapfrog_positions_3d + + subroutine warpx_push_leapfrog_positions_2d(particles, ns, np, & + vx_p, vy_p, dt, & + prob_lo, prob_hi) & + bind(c,name='warpx_push_leapfrog_positions_2d') + integer, value, intent(in) :: ns, np + real(amrex_real), intent(inout) :: particles(ns,np) + real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) + real(amrex_real), intent(in) :: dt + real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) + + integer n + + do n = 1, np + + particles(1, n) = particles(1, n) + dt * vx_p(n) + particles(2, n) = particles(2, n) + dt * vy_p(n) + +! bounce off the walls in the x... + do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) + if (particles(1, n) .lt. prob_lo(1)) then + particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) + else + particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) + end if + vx_p(n) = -vx_p(n) + end do + +! ... y... + do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) + if (particles(2, n) .lt. prob_lo(2)) then + particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) + else + particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) + end if + vy_p(n) = -vy_p(n) + end do + + end do + + end subroutine warpx_push_leapfrog_positions_2d + +end module warpx_ES_push_particles diff --git a/Source/Python/Make.package b/Source/Python/Make.package new file mode 100644 index 000000000..71bd4ebe8 --- /dev/null +++ b/Source/Python/Make.package @@ -0,0 +1,11 @@ +ifeq ($(USE_PYTHON_MAIN),TRUE) + CEXE_sources += WarpXWrappers.cpp + CEXE_headers += WarpXWrappers.h +endif +CEXE_sources += WarpXWrappers.cpp +CEXE_sources += WarpX_py.cpp +CEXE_headers += WarpXWrappers.h +CEXE_headers += WarpX_py.H + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Python +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Python diff --git a/Source/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 6ee06a7d5..6ee06a7d5 100644 --- a/Source/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp diff --git a/Source/WarpXWrappers.h b/Source/Python/WarpXWrappers.h index 07d6f80f7..07d6f80f7 100644 --- a/Source/WarpXWrappers.h +++ b/Source/Python/WarpXWrappers.h diff --git a/Source/WarpX_py.H b/Source/Python/WarpX_py.H index d8cf22155..d8cf22155 100644 --- a/Source/WarpX_py.H +++ b/Source/Python/WarpX_py.H diff --git a/Source/WarpX_py.cpp b/Source/Python/WarpX_py.cpp index 276d637d7..276d637d7 100644 --- a/Source/WarpX_py.cpp +++ b/Source/Python/WarpX_py.cpp diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package new file mode 100644 index 000000000..8cd124ea3 --- /dev/null +++ b/Source/Utils/Make.package @@ -0,0 +1,10 @@ +F90EXE_sources += utils_ES.F90 +CEXE_sources += WarpXMovingWindow.cpp +CEXE_sources += WarpXConst.cpp +CEXE_sources += WarpXTagging.cpp +CEXE_sources += WarpXUtil.cpp +CEXE_headers += WarpXConst.H +CEXE_headers += WarpXUtil.H + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/WarpXConst.H b/Source/Utils/WarpXConst.H index fdb92f666..fdb92f666 100644 --- a/Source/WarpXConst.H +++ b/Source/Utils/WarpXConst.H diff --git a/Source/WarpXConst.cpp b/Source/Utils/WarpXConst.cpp index 4ca729872..bd3ebc3ef 100644 --- a/Source/WarpXConst.cpp +++ b/Source/Utils/WarpXConst.cpp @@ -8,7 +8,7 @@ #include <WarpX.H> #include <WarpXConst.H> #include <WarpX_f.H> -#include <ParticleContainer.H> +#include <MultiParticleContainer.H> std::string UserConstants::replaceStringValue(std::string math_expr){ std::string pattern, value_str; diff --git a/Source/WarpXMove.cpp b/Source/Utils/WarpXMovingWindow.cpp index 908c70573..908c70573 100644 --- a/Source/WarpXMove.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp diff --git a/Source/WarpXTagging.cpp b/Source/Utils/WarpXTagging.cpp index 21acbefdc..21acbefdc 100644 --- a/Source/WarpXTagging.cpp +++ b/Source/Utils/WarpXTagging.cpp diff --git a/Source/WarpXUtil.H b/Source/Utils/WarpXUtil.H index a973c2283..a973c2283 100644 --- a/Source/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H diff --git a/Source/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4a884330a..4a884330a 100644 --- a/Source/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp diff --git a/Source/Utils/utils_ES.F90 b/Source/Utils/utils_ES.F90 new file mode 100644 index 000000000..887a5ef15 --- /dev/null +++ b/Source/Utils/utils_ES.F90 @@ -0,0 +1,191 @@ +module warpx_ES_utils + + use iso_c_binding + use amrex_fort_module, only : amrex_real + + implicit none + +contains + + subroutine warpx_sum_fine_to_crse_nodal_3d (lo, hi, lrat, crse, clo, chi, fine, flo, fhi) & + bind(c, name="warpx_sum_fine_to_crse_nodal_3d") + + integer, intent(in) :: lo(3), hi(3) + integer, intent(in) :: clo(3), chi(3) + integer, intent(in) :: flo(3), fhi(3) + integer, intent(in) :: lrat(3) + real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3)) + real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3)) + + integer :: i, j, k, ii, jj, kk + + do k = lo(3), hi(3) + kk = k * lrat(3) + do j = lo(2), hi(2) + jj = j * lrat(2) + do i = lo(1), hi(1) + ii = i * lrat(1) + crse(i,j,k) = fine(ii,jj,kk) + & +! These six fine nodes are shared by two coarse nodes... + 0.5d0 * (fine(ii-1,jj,kk) + fine(ii+1,jj,kk) + & + fine(ii,jj-1,kk) + fine(ii,jj+1,kk) + & + fine(ii,jj,kk-1) + fine(ii,jj,kk+1)) + & +! ... these twelve are shared by four... + 0.25d0 * (fine(ii,jj-1,kk-1) + fine(ii,jj+1,kk-1) + & + fine(ii,jj-1,kk+1) + fine(ii,jj+1,kk+1) + & + fine(ii-1,jj,kk-1) + fine(ii+1,jj,kk-1) + & + fine(ii-1,jj,kk+1) + fine(ii+1,jj,kk+1) + & + fine(ii-1,jj-1,kk) + fine(ii+1,jj-1,kk) + & + fine(ii-1,jj+1,kk) + fine(ii+1,jj+1,kk)) + & +! ... and these eight are shared by eight + 0.125d0 * (fine(ii-1,jj-1,kk-1) + fine(ii-1,jj-1,kk+1) + & + fine(ii-1,jj+1,kk-1) + fine(ii-1,jj+1,kk+1) + & + fine(ii+1,jj-1,kk-1) + fine(ii+1,jj-1,kk+1) + & + fine(ii+1,jj+1,kk-1) + fine(ii+1,jj+1,kk+1)) +! ... note that we have 27 nodes in total... + crse(i,j,k) = crse(i,j,k) / 8.d0 + end do + end do + end do + + end subroutine warpx_sum_fine_to_crse_nodal_3d + + subroutine warpx_sum_fine_to_crse_nodal_2d (lo, hi, lrat, crse, clo, chi, fine, flo, fhi) & + bind(c, name="warpx_sum_fine_to_crse_nodal_2d") + + integer, intent(in) :: lo(2), hi(2) + integer, intent(in) :: clo(2), chi(2) + integer, intent(in) :: flo(2), fhi(2) + integer, intent(in) :: lrat(2) + real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2)) + real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2)) + + integer :: i, j, ii, jj + + do j = lo(2), hi(2) + jj = j * lrat(2) + do i = lo(1), hi(1) + ii = i * lrat(1) + crse(i,j) = fine(ii,jj) + & +! These four fine nodes are shared by two coarse nodes... + 0.5d0 * (fine(ii-1,jj) + fine(ii+1,jj) + & + fine(ii,jj-1) + fine(ii,jj+1)) + & +! ... and these four are shared by four... + 0.25d0 * (fine(ii-1,jj-1) + fine(ii-1,jj+1) + & + fine(ii-1,jj+1) + fine(ii+1,jj+1)) +! ... note that we have 9 nodes in total... + crse(i,j) = crse(i,j) / 4.d0 + end do + end do + + end subroutine warpx_sum_fine_to_crse_nodal_2d + + subroutine warpx_zero_out_bndry_3d (lo, hi, input_data, bndry_data, mask) & + bind(c,name='warpx_zero_out_bndry_3d') + + integer(c_int), intent(in ) :: lo(3), hi(3) + double precision, intent(inout) :: input_data(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) + double precision, intent(inout) :: bndry_data(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) + integer(c_int), intent(in ) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) + + integer :: i, j, k + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + if (mask(i,j,k) .eq. 1) then + bndry_data(i,j,k) = input_data(i,j,k) + input_data(i,j,k) = 0.d0 + end if + end do + end do + end do + + end subroutine warpx_zero_out_bndry_3d + + subroutine warpx_zero_out_bndry_2d (lo, hi, input_data, bndry_data, mask) & + bind(c,name='warpx_zero_out_bndry_2d') + + integer(c_int), intent(in ) :: lo(2), hi(2) + double precision, intent(inout) :: input_data(lo(1):hi(1),lo(2):hi(2)) + double precision, intent(inout) :: bndry_data(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) + integer(c_int), intent(in ) :: mask (lo(1):hi(1),lo(2):hi(2)) + + integer :: i, j + + do j = lo(2), hi(2) + do i = lo(1), hi(1) + if (mask(i,j) .eq. 1) then + bndry_data(i,j) = input_data(i,j) + input_data(i,j) = 0.d0 + end if + end do + end do + + end subroutine warpx_zero_out_bndry_2d + + subroutine warpx_build_mask_3d (lo, hi, tmp_mask, mask, ncells) & + bind(c,name='warpx_build_mask_3d') + integer(c_int), intent(in ) :: lo(3), hi(3) + integer(c_int), intent(in ) :: ncells + integer(c_int), intent(in ) :: tmp_mask(lo(1)-ncells:hi(1)+ncells,lo(2)-ncells:hi(2)+ncells,lo(3)-ncells:hi(3)+ncells) + integer(c_int), intent(inout) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) + + integer :: i, j, k, ii, jj, kk, total + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + total = 0 + do ii = i-ncells, i+ncells + do jj = j-ncells, j+ncells + do kk = k-ncells, k+ncells + total = total + tmp_mask(ii, jj, kk) + end do + end do + end do + + if (total .gt. 0) then + mask(i,j,k) = 1 + else + mask(i,j,k) = 0 + end if + + end do + end do + end do + + end subroutine warpx_build_mask_3d + + subroutine warpx_build_mask_2d (lo, hi, tmp_mask, mask, ncells) & + bind(c,name='warpx_build_mask_2d') + integer(c_int), intent(in ) :: lo(2), hi(2) + integer(c_int), intent(in ) :: ncells + integer(c_int), intent(in ) :: tmp_mask(lo(1)-ncells:hi(1)+ncells,lo(2)-ncells:hi(2)+ncells) + integer(c_int), intent(inout) :: mask (lo(1):hi(1),lo(2):hi(2)) + + integer :: i, j, ii, jj, total + + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + total = 0 + do ii = i-ncells, i+ncells + do jj = j-ncells, j+ncells + total = total + tmp_mask(ii, jj) + end do + end do + + if (total .gt. 0) then + mask(i,j) = 1 + else + mask(i,j) = 0 + end if + + end do + end do + + end subroutine warpx_build_mask_2d + +end module warpx_ES_utils diff --git a/Source/WarpX.H b/Source/WarpX.H index bf82bceca..2462182ed 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -19,9 +19,9 @@ #include <AMReX_Interpolater.H> #include <AMReX_FillPatchUtil.H> -#include <ParticleContainer.H> -#include <WarpXPML.H> -#include <WarpXBoostedFrameDiagnostic.H> +#include <MultiParticleContainer.H> +#include <PML.H> +#include <BoostedFrameDiagnostic.H> #ifdef WARPX_USE_PSATD #include <fftw3.h> @@ -501,7 +501,7 @@ private: bool plot_part_per_proc = false; bool plot_proc_number = false; bool plot_dive = false; - bool plot_divb = true; + bool plot_divb = false; bool plot_rho = false; bool plot_F = false; bool plot_finepatch = false; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index ba54acee0..bf0d9b255 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -851,6 +851,21 @@ WarpX::RefRatio (int lev) } void +WarpX::Evolve (int numsteps) { + BL_PROFILE_REGION("WarpX::Evolve()"); + +#ifdef WARPX_DO_ELECTROSTATIC + if (do_electrostatic) { + EvolveES(numsteps); + } else { + EvolveEM(numsteps); + } +#else + EvolveEM(numsteps); +#endif // WARPX_DO_ELECTROSTATIC +} + +void WarpX::ComputeDivB (MultiFab& divB, int dcomp, const std::array<const MultiFab*, 3>& B, const std::array<Real,3>& dx) diff --git a/Source/WarpXProb.cpp b/Source/WarpXProb.cpp deleted file mode 100644 index ecbf1afaf..000000000 --- a/Source/WarpXProb.cpp +++ /dev/null @@ -1,78 +0,0 @@ -#include <WarpX.H> -#include <AMReX_ParmParse.H> - -using namespace amrex; - -void -WarpX::InitLevelData (int lev, Real time) -{ - for (int i = 0; i < 3; ++i) { - current_fp[lev][i]->setVal(0.0); - Efield_fp[lev][i]->setVal(0.0); - Bfield_fp[lev][i]->setVal(0.0); - } - - if (lev > 0) { - for (int i = 0; i < 3; ++i) { - Efield_aux[lev][i]->setVal(0.0); - Bfield_aux[lev][i]->setVal(0.0); - - current_cp[lev][i]->setVal(0.0); - Efield_cp[lev][i]->setVal(0.0); - Bfield_cp[lev][i]->setVal(0.0); - } - } - - if (F_fp[lev]) { - F_fp[lev]->setVal(0.0); - } - - if (rho_fp[lev]) { - rho_fp[lev]->setVal(0.0); - } - - if (F_cp[lev]) { - F_cp[lev]->setVal(0.0); - } - - if (rho_cp[lev]) { - rho_cp[lev]->setVal(0.0); - } - - if (costs[lev]) { - costs[lev]->setVal(0.0); - } -} - -#ifdef WARPX_USE_PSATD - -void -WarpX::InitLevelDataFFT (int lev, Real time) -{ - Efield_fp_fft[lev][0]->setVal(0.0); - Efield_fp_fft[lev][1]->setVal(0.0); - Efield_fp_fft[lev][2]->setVal(0.0); - Bfield_fp_fft[lev][0]->setVal(0.0); - Bfield_fp_fft[lev][1]->setVal(0.0); - Bfield_fp_fft[lev][2]->setVal(0.0); - current_fp_fft[lev][0]->setVal(0.0); - current_fp_fft[lev][1]->setVal(0.0); - current_fp_fft[lev][2]->setVal(0.0); - rho_fp_fft[lev]->setVal(0.0); - - if (lev > 0) - { - Efield_cp_fft[lev][0]->setVal(0.0); - Efield_cp_fft[lev][1]->setVal(0.0); - Efield_cp_fft[lev][2]->setVal(0.0); - Bfield_cp_fft[lev][0]->setVal(0.0); - Bfield_cp_fft[lev][1]->setVal(0.0); - Bfield_cp_fft[lev][2]->setVal(0.0); - current_cp_fft[lev][0]->setVal(0.0); - current_cp_fft[lev][1]->setVal(0.0); - current_cp_fft[lev][2]->setVal(0.0); - rho_cp_fft[lev]->setVal(0.0); - } -} - -#endif diff --git a/Source/WarpX_electrostatic.F90 b/Source/WarpX_electrostatic.F90 deleted file mode 100644 index 1e5140765..000000000 --- a/Source/WarpX_electrostatic.F90 +++ /dev/null @@ -1,990 +0,0 @@ -module warpx_electrostatic_module - - use iso_c_binding - use amrex_fort_module, only : amrex_real - - implicit none - -contains - - subroutine warpx_sum_fine_to_crse_nodal_3d (lo, hi, lrat, crse, clo, chi, fine, flo, fhi) & - bind(c, name="warpx_sum_fine_to_crse_nodal_3d") - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: clo(3), chi(3) - integer, intent(in) :: flo(3), fhi(3) - integer, intent(in) :: lrat(3) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3)) - real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3)) - - integer :: i, j, k, ii, jj, kk - - do k = lo(3), hi(3) - kk = k * lrat(3) - do j = lo(2), hi(2) - jj = j * lrat(2) - do i = lo(1), hi(1) - ii = i * lrat(1) - crse(i,j,k) = fine(ii,jj,kk) + & -! These six fine nodes are shared by two coarse nodes... - 0.5d0 * (fine(ii-1,jj,kk) + fine(ii+1,jj,kk) + & - fine(ii,jj-1,kk) + fine(ii,jj+1,kk) + & - fine(ii,jj,kk-1) + fine(ii,jj,kk+1)) + & -! ... these twelve are shared by four... - 0.25d0 * (fine(ii,jj-1,kk-1) + fine(ii,jj+1,kk-1) + & - fine(ii,jj-1,kk+1) + fine(ii,jj+1,kk+1) + & - fine(ii-1,jj,kk-1) + fine(ii+1,jj,kk-1) + & - fine(ii-1,jj,kk+1) + fine(ii+1,jj,kk+1) + & - fine(ii-1,jj-1,kk) + fine(ii+1,jj-1,kk) + & - fine(ii-1,jj+1,kk) + fine(ii+1,jj+1,kk)) + & -! ... and these eight are shared by eight - 0.125d0 * (fine(ii-1,jj-1,kk-1) + fine(ii-1,jj-1,kk+1) + & - fine(ii-1,jj+1,kk-1) + fine(ii-1,jj+1,kk+1) + & - fine(ii+1,jj-1,kk-1) + fine(ii+1,jj-1,kk+1) + & - fine(ii+1,jj+1,kk-1) + fine(ii+1,jj+1,kk+1)) -! ... note that we have 27 nodes in total... - crse(i,j,k) = crse(i,j,k) / 8.d0 - end do - end do - end do - - end subroutine warpx_sum_fine_to_crse_nodal_3d - - subroutine warpx_sum_fine_to_crse_nodal_2d (lo, hi, lrat, crse, clo, chi, fine, flo, fhi) & - bind(c, name="warpx_sum_fine_to_crse_nodal_2d") - - integer, intent(in) :: lo(2), hi(2) - integer, intent(in) :: clo(2), chi(2) - integer, intent(in) :: flo(2), fhi(2) - integer, intent(in) :: lrat(2) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2)) - real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2)) - - integer :: i, j, ii, jj - - do j = lo(2), hi(2) - jj = j * lrat(2) - do i = lo(1), hi(1) - ii = i * lrat(1) - crse(i,j) = fine(ii,jj) + & -! These four fine nodes are shared by two coarse nodes... - 0.5d0 * (fine(ii-1,jj) + fine(ii+1,jj) + & - fine(ii,jj-1) + fine(ii,jj+1)) + & -! ... and these four are shared by four... - 0.25d0 * (fine(ii-1,jj-1) + fine(ii-1,jj+1) + & - fine(ii-1,jj+1) + fine(ii+1,jj+1)) -! ... note that we have 9 nodes in total... - crse(i,j) = crse(i,j) / 4.d0 - end do - end do - - end subroutine warpx_sum_fine_to_crse_nodal_2d - - subroutine warpx_zero_out_bndry_3d (lo, hi, input_data, bndry_data, mask) & - bind(c,name='warpx_zero_out_bndry_3d') - - integer(c_int), intent(in ) :: lo(3), hi(3) - double precision, intent(inout) :: input_data(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) - double precision, intent(inout) :: bndry_data(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) - integer(c_int), intent(in ) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) - - integer :: i, j, k - - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - if (mask(i,j,k) .eq. 1) then - bndry_data(i,j,k) = input_data(i,j,k) - input_data(i,j,k) = 0.d0 - end if - end do - end do - end do - - end subroutine warpx_zero_out_bndry_3d - - subroutine warpx_zero_out_bndry_2d (lo, hi, input_data, bndry_data, mask) & - bind(c,name='warpx_zero_out_bndry_2d') - - integer(c_int), intent(in ) :: lo(2), hi(2) - double precision, intent(inout) :: input_data(lo(1):hi(1),lo(2):hi(2)) - double precision, intent(inout) :: bndry_data(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) - integer(c_int), intent(in ) :: mask (lo(1):hi(1),lo(2):hi(2)) - - integer :: i, j - - do j = lo(2), hi(2) - do i = lo(1), hi(1) - if (mask(i,j) .eq. 1) then - bndry_data(i,j) = input_data(i,j) - input_data(i,j) = 0.d0 - end if - end do - end do - - end subroutine warpx_zero_out_bndry_2d - - subroutine warpx_build_mask_3d (lo, hi, tmp_mask, mask, ncells) & - bind(c,name='warpx_build_mask_3d') - integer(c_int), intent(in ) :: lo(3), hi(3) - integer(c_int), intent(in ) :: ncells - integer(c_int), intent(in ) :: tmp_mask(lo(1)-ncells:hi(1)+ncells,lo(2)-ncells:hi(2)+ncells,lo(3)-ncells:hi(3)+ncells) - integer(c_int), intent(inout) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) - - integer :: i, j, k, ii, jj, kk, total - - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - - total = 0 - do ii = i-ncells, i+ncells - do jj = j-ncells, j+ncells - do kk = k-ncells, k+ncells - total = total + tmp_mask(ii, jj, kk) - end do - end do - end do - - if (total .gt. 0) then - mask(i,j,k) = 1 - else - mask(i,j,k) = 0 - end if - - end do - end do - end do - - end subroutine warpx_build_mask_3d - - subroutine warpx_build_mask_2d (lo, hi, tmp_mask, mask, ncells) & - bind(c,name='warpx_build_mask_2d') - integer(c_int), intent(in ) :: lo(2), hi(2) - integer(c_int), intent(in ) :: ncells - integer(c_int), intent(in ) :: tmp_mask(lo(1)-ncells:hi(1)+ncells,lo(2)-ncells:hi(2)+ncells) - integer(c_int), intent(inout) :: mask (lo(1):hi(1),lo(2):hi(2)) - - integer :: i, j, ii, jj, total - - do j = lo(2), hi(2) - do i = lo(1), hi(1) - - total = 0 - do ii = i-ncells, i+ncells - do jj = j-ncells, j+ncells - total = total + tmp_mask(ii, jj) - end do - end do - - if (total .gt. 0) then - mask(i,j) = 1 - else - mask(i,j) = 0 - end if - - end do - end do - - end subroutine warpx_build_mask_2d - -! This routine computes the node-centered electric field given a node-centered phi. -! The gradient is computed using 2nd-order centered differences. It assumes the -! Boundary conditions have already been set and that you have two rows of ghost cells -! for phi and one row of ghost cells for Ex, Ey, and Ez. -! Note that this routine includes the minus sign in E = - grad phi. -! -! Arguments: -! lo, hi: The corners of the valid box over which the gradient is taken -! Ex, Ey, Ez: The electric field in the x, y, and z directions. -! dx: The cell spacing -! - subroutine warpx_compute_E_nodal_3d (lo, hi, phi, Ex, Ey, Ez, dx) & - bind(c,name='warpx_compute_E_nodal_3d') - integer(c_int), intent(in) :: lo(3), hi(3) - real(amrex_real), intent(in) :: dx(3) - real(amrex_real), intent(in ) :: phi(lo(1)-2:hi(1)+2,lo(2)-2:hi(2)+2,lo(3)-2:hi(3)+2) - real(amrex_real), intent(inout) :: Ex (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) - real(amrex_real), intent(inout) :: Ey (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) - real(amrex_real), intent(inout) :: Ez (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) - - integer :: i, j, k - real(amrex_real) :: fac(3) - - fac = 0.5d0 / dx - - do k = lo(3)-1, hi(3)+1 - do j = lo(2)-1, hi(2)+1 - do i = lo(1)-1, hi(1)+1 - - Ex(i,j,k) = fac(1) * (phi(i-1,j,k) - phi(i+1,j,k)) - Ey(i,j,k) = fac(2) * (phi(i,j-1,k) - phi(i,j+1,k)) - Ez(i,j,k) = fac(3) * (phi(i,j,k-1) - phi(i,j,k+1)) - - end do - end do - end do - - end subroutine warpx_compute_E_nodal_3d - - subroutine warpx_compute_E_nodal_2d (lo, hi, phi, Ex, Ey, dx) & - bind(c,name='warpx_compute_E_nodal_2d') - integer(c_int), intent(in) :: lo(2), hi(2) - real(amrex_real), intent(in) :: dx(2) - real(amrex_real), intent(in ) :: phi(lo(1)-2:hi(1)+2,lo(2)-2:hi(2)+2) - real(amrex_real), intent(inout) :: Ex (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) - real(amrex_real), intent(inout) :: Ey (lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) - - integer :: i, j - real(amrex_real) :: fac(2) - - fac = 0.5d0 / dx - - do j = lo(2)-1, hi(2)+1 - do i = lo(1)-1, hi(1)+1 - - Ex(i,j) = fac(1) * (phi(i-1,j) - phi(i+1,j)) - Ey(i,j) = fac(2) * (phi(i,j-1) - phi(i,j+1)) - - end do - end do - - end subroutine warpx_compute_E_nodal_2d - - -! This routine computes the charge density due to the particles using cloud-in-cell -! deposition. The Fab rho is assumed to be node-centered. -! -! Arguments: -! particles : a pointer to the particle array-of-structs -! ns : the stride length of particle struct (the size of the struct in number of reals) -! np : the number of particles -! weights : the particle weights -! charge : the charge of this particle species -! rho : a Fab that will contain the charge density on exit -! lo : a pointer to the lo corner of this valid box for rho, in index space -! hi : a pointer to the hi corner of this valid box for rho, in index space -! plo : the real position of the left-hand corner of the problem domain -! dx : the mesh spacing -! ng : the number of ghost cells in rho -! - subroutine warpx_deposit_cic_3d(particles, ns, np, & - weights, charge, rho, lo, hi, plo, dx, & - ng) & - bind(c,name='warpx_deposit_cic_3d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(in) :: weights(np) - real(amrex_real), intent(in) :: charge - integer, intent(in) :: lo(3) - integer, intent(in) :: hi(3) - integer, intent(in) :: ng - real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: plo(3) - real(amrex_real), intent(in) :: dx(3) - - integer i, j, k, n - real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi - real(amrex_real) lx, ly, lz - real(amrex_real) inv_dx(3) - real(amrex_real) qp, inv_vol - - inv_dx = 1.0d0/dx - - inv_vol = inv_dx(1) * inv_dx(2) * inv_dx(3) - - do n = 1, np - - qp = weights(n) * charge * inv_vol - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - lz = (particles(3, n) - plo(3))*inv_dx(3) - - i = floor(lx) - j = floor(ly) - k = floor(lz) - - wx_hi = lx - i - wy_hi = ly - j - wz_hi = lz - k - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - wz_lo = 1.0d0 - wz_hi - - rho(i, j, k ) = rho(i, j, k ) + wx_lo*wy_lo*wz_lo*qp - rho(i, j, k+1) = rho(i, j, k+1) + wx_lo*wy_lo*wz_hi*qp - rho(i, j+1, k ) = rho(i, j+1, k ) + wx_lo*wy_hi*wz_lo*qp - rho(i, j+1, k+1) = rho(i, j+1, k+1) + wx_lo*wy_hi*wz_hi*qp - rho(i+1, j, k ) = rho(i+1, j, k ) + wx_hi*wy_lo*wz_lo*qp - rho(i+1, j, k+1) = rho(i+1, j, k+1) + wx_hi*wy_lo*wz_hi*qp - rho(i+1, j+1, k ) = rho(i+1, j+1, k ) + wx_hi*wy_hi*wz_lo*qp - rho(i+1, j+1, k+1) = rho(i+1, j+1, k+1) + wx_hi*wy_hi*wz_hi*qp - - end do - - end subroutine warpx_deposit_cic_3d - - subroutine warpx_deposit_cic_2d(particles, ns, np, & - weights, charge, rho, lo, hi, plo, dx, & - ng) & - bind(c,name='warpx_deposit_cic_2d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(in) :: weights(np) - real(amrex_real), intent(in) :: charge - integer, intent(in) :: lo(2) - integer, intent(in) :: hi(2) - integer, intent(in) :: ng - real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: plo(2) - real(amrex_real), intent(in) :: dx(2) - - integer i, j, n - real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi - real(amrex_real) lx, ly - real(amrex_real) inv_dx(2) - real(amrex_real) qp, inv_vol - - inv_dx = 1.0d0/dx - - inv_vol = inv_dx(1) * inv_dx(2) - - do n = 1, np - - qp = weights(n) * charge * inv_vol - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - - i = floor(lx) - j = floor(ly) - - wx_hi = lx - i - wy_hi = ly - j - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - - rho(i, j ) = rho(i, j ) + wx_lo*wy_lo*qp - rho(i, j+1) = rho(i, j+1) + wx_lo*wy_hi*qp - rho(i+1, j ) = rho(i+1, j ) + wx_hi*wy_lo*qp - rho(i+1, j+1) = rho(i+1, j+1) + wx_hi*wy_hi*qp - - end do - - end subroutine warpx_deposit_cic_2d - - -! This routine interpolates the electric field to the particle positions -! using cloud-in-cell interpolation. The electric fields are assumed to be -! node-centered. -! -! Arguments: -! particles : a pointer to the particle array-of-structs -! ns : the stride length of particle struct (the size of the struct in number of reals) -! np : the number of particles -! Ex_p : the electric field in the x-direction at the particle positions (output) -! Ey_p : the electric field in the y-direction at the particle positions (output) -! Ez_p : the electric field in the z-direction at the particle positions (output) -! Ex, Ey, Ez: Fabs conting the electric field on the mesh -! lo : a pointer to the lo corner of this valid box, in index space -! hi : a pointer to the hi corner of this valid box, in index space -! plo : the real position of the left-hand corner of the problem domain -! dx : the mesh spacing -! ng : the number of ghost cells for the E-field -! - subroutine warpx_interpolate_cic_3d(particles, ns, np, & - Ex_p, Ey_p, Ez_p, & - Ex, Ey, Ez, & - lo, hi, plo, dx, ng) & - bind(c,name='warpx_interpolate_cic_3d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) - integer, intent(in) :: ng - integer, intent(in) :: lo(3) - integer, intent(in) :: hi(3) - real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: Ez(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: plo(3) - real(amrex_real), intent(in) :: dx(3) - - integer i, j, k, n - real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi - real(amrex_real) lx, ly, lz - real(amrex_real) inv_dx(3) - inv_dx = 1.0d0/dx - - do n = 1, np - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - lz = (particles(3, n) - plo(3))*inv_dx(3) - - i = floor(lx) - j = floor(ly) - k = floor(lz) - - wx_hi = lx - i - wy_hi = ly - j - wz_hi = lz - k - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - wz_lo = 1.0d0 - wz_hi - - Ex_p(n) = wx_lo*wy_lo*wz_lo*Ex(i, j, k ) + & - wx_lo*wy_lo*wz_hi*Ex(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*Ex(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*Ex(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*Ex(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) - - Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & - wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*Ey(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*Ey(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) - - Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & - wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*Ez(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*Ez(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) - - end do - - end subroutine warpx_interpolate_cic_3d - - - subroutine warpx_interpolate_cic_2d(particles, ns, np, & - Ex_p, Ey_p, & - Ex, Ey, & - lo, hi, plo, dx, ng) & - bind(c,name='warpx_interpolate_cic_2d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) - integer, intent(in) :: ng - integer, intent(in) :: lo(2) - integer, intent(in) :: hi(2) - real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: plo(2) - real(amrex_real), intent(in) :: dx(2) - - integer i, j, n - real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi - real(amrex_real) lx, ly - real(amrex_real) inv_dx(2) - inv_dx = 1.0d0/dx - - do n = 1, np - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - - i = floor(lx) - j = floor(ly) - - wx_hi = lx - i - wy_hi = ly - j - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - - Ex_p(n) = wx_lo*wy_lo*Ex(i, j ) + & - wx_lo*wy_hi*Ex(i, j+1) + & - wx_hi*wy_lo*Ex(i+1, j ) + & - wx_hi*wy_hi*Ex(i+1, j+1) - - Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & - wx_lo*wy_hi*Ey(i, j+1) + & - wx_hi*wy_lo*Ey(i+1, j ) + & - wx_hi*wy_hi*Ey(i+1, j+1) - - end do - - end subroutine warpx_interpolate_cic_2d - - - subroutine warpx_interpolate_cic_two_levels_3d(particles, ns, np, & - Ex_p, Ey_p, Ez_p, & - Ex, Ey, Ez, & - lo, hi, dx, & - cEx, cEy, cEz, & - mask, & - clo, chi, cdx, & - plo, ng, lev) & - bind(c,name='warpx_interpolate_cic_two_levels_3d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) - integer, intent(in) :: ng, lev - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: clo(3), chi(3) - real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: Ez(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: cEx(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) - real(amrex_real), intent(in) :: cEy(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) - real(amrex_real), intent(in) :: cEz(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng, clo(3)-ng:chi(3)+ng) - integer(c_int), intent(in) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) - real(amrex_real), intent(in) :: plo(3) - real(amrex_real), intent(in) :: dx(3), cdx(3) - - integer i, j, k, n - real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi - real(amrex_real) lx, ly, lz - real(amrex_real) inv_dx(3), inv_cdx(3) - inv_dx = 1.0d0/dx - inv_cdx = 1.0d0/cdx - - do n = 1, np - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - lz = (particles(3, n) - plo(3))*inv_dx(3) - - i = floor(lx) - j = floor(ly) - k = floor(lz) - -! use the coarse E if near the level boundary - if (lev .eq. 1 .and. mask(i,j,k) .eq. 1) then - - lx = (particles(1, n) - plo(1))*inv_cdx(1) - ly = (particles(2, n) - plo(2))*inv_cdx(2) - lz = (particles(3, n) - plo(3))*inv_cdx(3) - - i = floor(lx) - j = floor(ly) - k = floor(lz) - - wx_hi = lx - i - wy_hi = ly - j - wz_hi = lz - k - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - wz_lo = 1.0d0 - wz_hi - - Ex_p(n) = wx_lo*wy_lo*wz_lo*cEx(i, j, k ) + & - wx_lo*wy_lo*wz_hi*cEx(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*cEx(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*cEx(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*cEx(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*cEx(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*cEx(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*cEx(i+1, j+1, k+1) - - Ey_p(n) = wx_lo*wy_lo*wz_lo*cEy(i, j, k ) + & - wx_lo*wy_lo*wz_hi*cEy(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*cEy(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*cEy(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*cEy(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*cEy(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*cEy(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*cEy(i+1, j+1, k+1) - - Ez_p(n) = wx_lo*wy_lo*wz_lo*cEz(i, j, k ) + & - wx_lo*wy_lo*wz_hi*cEz(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*cEz(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*cEz(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*cEz(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*cEz(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*cEz(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*cEz(i+1, j+1, k+1) - -! otherwise use the fine - else - - wx_hi = lx - i - wy_hi = ly - j - wz_hi = lz - k - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - wz_lo = 1.0d0 - wz_hi - - Ex_p(n) = wx_lo*wy_lo*wz_lo*Ex(i, j, k ) + & - wx_lo*wy_lo*wz_hi*Ex(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*Ex(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*Ex(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*Ex(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) - - Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & - wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*Ey(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*Ey(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) - - Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & - wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & - wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & - wx_lo*wy_hi*wz_hi*Ez(i, j+1, k+1) + & - wx_hi*wy_lo*wz_lo*Ez(i+1, j, k ) + & - wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & - wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & - wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) - - end if - - end do - - end subroutine warpx_interpolate_cic_two_levels_3d - - - subroutine warpx_interpolate_cic_two_levels_2d(particles, ns, np, & - Ex_p, Ey_p, & - Ex, Ey, & - lo, hi, dx, & - cEx, cEy, & - mask, & - clo, chi, cdx, & - plo, ng, lev) & - bind(c,name='warpx_interpolate_cic_two_levels_2d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) - integer, intent(in) :: ng, lev - integer, intent(in) :: lo(2), hi(2) - integer, intent(in) :: clo(2), chi(2) - real(amrex_real), intent(in) :: Ex(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: Ey(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: cEx(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng) - real(amrex_real), intent(in) :: cEy(clo(1)-ng:chi(1)+ng, clo(2)-ng:chi(2)+ng) - integer(c_int), intent(in) :: mask (lo(1):hi(1),lo(2):hi(2)) - real(amrex_real), intent(in) :: plo(2) - real(amrex_real), intent(in) :: dx(2), cdx(2) - - integer i, j, n - real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi - real(amrex_real) lx, ly - real(amrex_real) inv_dx(2), inv_cdx(2) - inv_dx = 1.0d0/dx - inv_cdx = 1.0d0/cdx - - do n = 1, np - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - - i = floor(lx) - j = floor(ly) - -! use the coarse E if near the level boundary - if (lev .eq. 1 .and. mask(i,j) .eq. 1) then - - lx = (particles(1, n) - plo(1))*inv_cdx(1) - ly = (particles(2, n) - plo(2))*inv_cdx(2) - - i = floor(lx) - j = floor(ly) - - wx_hi = lx - i - wy_hi = ly - j - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - - Ex_p(n) = wx_lo*wy_lo*cEx(i, j ) + & - wx_lo*wy_hi*cEx(i, j+1) + & - wx_hi*wy_lo*cEx(i+1, j ) + & - wx_hi*wy_hi*cEx(i+1, j+1) - - Ey_p(n) = wx_lo*wy_lo*cEy(i, j ) + & - wx_lo*wy_hi*cEy(i, j+1) + & - wx_hi*wy_lo*cEy(i+1, j ) + & - wx_hi*wy_hi*cEy(i+1, j+1) - -! otherwise use the fine - else - - wx_hi = lx - i - wy_hi = ly - j - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - - Ex_p(n) = wx_lo*wy_lo*Ex(i, j ) + & - wx_lo*wy_hi*Ex(i, j+1) + & - wx_hi*wy_lo*Ex(i+1, j ) + & - wx_hi*wy_hi*Ex(i+1, j+1) - - Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & - wx_lo*wy_hi*Ey(i, j+1) + & - wx_hi*wy_lo*Ey(i+1, j ) + & - wx_hi*wy_hi*Ey(i+1, j+1) - - end if - - end do - - end subroutine warpx_interpolate_cic_two_levels_2d - - -! -! This routine updates the particle positions and velocities using the -! leapfrog time integration algorithm, given the electric fields at the -! particle positions. It also enforces specular reflection off the domain -! walls. -! -! Arguments: -! particles : a pointer to the particle array-of-structs -! ns : the stride length of particle struct (the size of the struct in number of reals) -! np : the number of particles -! vx_p : the particle x-velocities -! vy_p : the particle y-velocities -! vz_p : the particle z-velocities -! Ex_p : the electric field in the x-direction at the particle positions -! Ey_p : the electric field in the y-direction at the particle positions -! Ez_p : the electric field in the z-direction at the particle positions -! charge : the charge of this particle species -! mass : the mass of this particle species -! dt : the time step -! prob_lo : the left-hand corner of the problem domain -! prob_hi : the right-hand corner of the problem domain -! - subroutine warpx_push_leapfrog_3d(particles, ns, np, & - vx_p, vy_p, vz_p, & - Ex_p, Ey_p, Ez_p, & - charge, mass, dt, & - prob_lo, prob_hi) & - bind(c,name='warpx_push_leapfrog_3d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) - real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np), Ez_p(np) - real(amrex_real), intent(in) :: charge - real(amrex_real), intent(in) :: mass - real(amrex_real), intent(in) :: dt - real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) - - integer n - real(amrex_real) fac - - fac = charge * dt / mass - - do n = 1, np - - vx_p(n) = vx_p(n) + fac * Ex_p(n) - vy_p(n) = vy_p(n) + fac * Ey_p(n) - vz_p(n) = vz_p(n) + fac * Ez_p(n) - - particles(1, n) = particles(1, n) + dt * vx_p(n) - particles(2, n) = particles(2, n) + dt * vy_p(n) - particles(3, n) = particles(3, n) + dt * vz_p(n) - -! bounce off the walls in the x... - do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) - if (particles(1, n) .lt. prob_lo(1)) then - particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) - else - particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) - end if - vx_p(n) = -vx_p(n) - end do - -! ... y... - do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) - if (particles(2, n) .lt. prob_lo(2)) then - particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) - else - particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) - end if - vy_p(n) = -vy_p(n) - end do - -! ... and z directions - do while (particles(3, n) .lt. prob_lo(3) .or. particles(3, n) .gt. prob_hi(3)) - if (particles(3, n) .lt. prob_lo(3)) then - particles(3, n) = 2.d0*prob_lo(3) - particles(3, n) - else - particles(3, n) = 2.d0*prob_hi(3) - particles(3, n) - end if - vz_p(n) = -vz_p(n) - end do - - end do - - end subroutine warpx_push_leapfrog_3d - - subroutine warpx_push_leapfrog_2d(particles, ns, np, & - vx_p, vy_p, & - Ex_p, Ey_p, & - charge, mass, dt, & - prob_lo, prob_hi) & - bind(c,name='warpx_push_leapfrog_2d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) - real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np) - real(amrex_real), intent(in) :: charge - real(amrex_real), intent(in) :: mass - real(amrex_real), intent(in) :: dt - real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) - - integer n - real(amrex_real) fac - - fac = charge * dt / mass - - do n = 1, np - - vx_p(n) = vx_p(n) + fac * Ex_p(n) - vy_p(n) = vy_p(n) + fac * Ey_p(n) - - particles(1, n) = particles(1, n) + dt * vx_p(n) - particles(2, n) = particles(2, n) + dt * vy_p(n) - -! bounce off the walls in the x... - do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) - if (particles(1, n) .lt. prob_lo(1)) then - particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) - else - particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) - end if - vx_p(n) = -vx_p(n) - end do - -! ... y... - do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) - if (particles(2, n) .lt. prob_lo(2)) then - particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) - else - particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) - end if - vy_p(n) = -vy_p(n) - end do - - end do - - end subroutine warpx_push_leapfrog_2d - - -! -! This routine advances the particle positions using the current -! velocity. This is needed to desynchronize the particle positions -! from the velocities after particle initialization. -! -! Arguments: -! particles : a pointer to the particle array-of-structs -! ns : the stride length of particle struct (the size of the struct in number of reals) -! np : the number of particles -! xx_p : the electric field in the x-direction at the particle positions -! vy_p : the electric field in the y-direction at the particle positions -! vz_p : the electric field in the z-direction at the particle positions -! dt : the time step -! prob_lo : the left-hand corner of the problem domain -! prob_hi : the right-hand corner of the problem domain -! - subroutine warpx_push_leapfrog_positions_3d(particles, ns, np, & - vx_p, vy_p, vz_p, dt, & - prob_lo, prob_hi) & - bind(c,name='warpx_push_leapfrog_positions_3d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) - real(amrex_real), intent(in) :: dt - real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) - - integer n - - do n = 1, np - - particles(1, n) = particles(1, n) + dt * vx_p(n) - particles(2, n) = particles(2, n) + dt * vy_p(n) - particles(3, n) = particles(3, n) + dt * vz_p(n) - -! bounce off the walls in the x... - do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) - if (particles(1, n) .lt. prob_lo(1)) then - particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) - else - particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) - end if - vx_p(n) = -vx_p(n) - end do - -! ... y... - do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) - if (particles(2, n) .lt. prob_lo(2)) then - particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) - else - particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) - end if - vy_p(n) = -vy_p(n) - end do - -! ... and z directions - do while (particles(3, n) .lt. prob_lo(3) .or. particles(3, n) .gt. prob_hi(3)) - if (particles(3, n) .lt. prob_lo(3)) then - particles(3, n) = 2.d0*prob_lo(3) - particles(3, n) - else - particles(3, n) = 2.d0*prob_hi(3) - particles(3, n) - end if - vz_p(n) = -vz_p(n) - end do - - end do - - end subroutine warpx_push_leapfrog_positions_3d - - subroutine warpx_push_leapfrog_positions_2d(particles, ns, np, & - vx_p, vy_p, dt, & - prob_lo, prob_hi) & - bind(c,name='warpx_push_leapfrog_positions_2d') - integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) - real(amrex_real), intent(in) :: dt - real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) - - integer n - - do n = 1, np - - particles(1, n) = particles(1, n) + dt * vx_p(n) - particles(2, n) = particles(2, n) + dt * vy_p(n) - -! bounce off the walls in the x... - do while (particles(1, n) .lt. prob_lo(1) .or. particles(1, n) .gt. prob_hi(1)) - if (particles(1, n) .lt. prob_lo(1)) then - particles(1, n) = 2.d0*prob_lo(1) - particles(1, n) - else - particles(1, n) = 2.d0*prob_hi(1) - particles(1, n) - end if - vx_p(n) = -vx_p(n) - end do - -! ... y... - do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) - if (particles(2, n) .lt. prob_lo(2)) then - particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) - else - particles(2, n) = 2.d0*prob_hi(2) - particles(2, n) - end if - vy_p(n) = -vy_p(n) - end do - - end do - - end subroutine warpx_push_leapfrog_positions_2d - -end module warpx_electrostatic_module diff --git a/tests/CurrentDeposition/GNUmakefile b/tests/CurrentDeposition/GNUmakefile deleted file mode 100644 index 05ac27a4b..000000000 --- a/tests/CurrentDeposition/GNUmakefile +++ /dev/null @@ -1,25 +0,0 @@ -AMREX_HOME ?= ../../../amrex -PICSAR_HOME ?= ../../../picsar - -USE_PARTICLES = TRUE - -DEBUG = FALSE -USE_MPI = TRUE -USE_OMP = TRUE -PROFILE = FALSE -COMP = gnu -DIM = 3 -PRECISION = DOUBLE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package -include $(AMREX_HOME)/Src/Base/Make.package -include $(PICSAR_HOME)/src/Make.package - -DEFINES += -DWARPX - -default: $(executable) - @echo SUCCESS - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/tests/CurrentDeposition/Make.package b/tests/CurrentDeposition/Make.package deleted file mode 100644 index 3129dcc37..000000000 --- a/tests/CurrentDeposition/Make.package +++ /dev/null @@ -1,9 +0,0 @@ - -CEXE_sources += main.cpp - -CEXE_headers += WarpX_f.H WarpXConst.H - -F90EXE_sources += WarpX_picsar.F90 - -INCLUDE_LOCATIONS += ../../Source -VPATH_LOCATIONS += ../../Source diff --git a/tests/CurrentDeposition/inputs b/tests/CurrentDeposition/inputs deleted file mode 100644 index 0be4f2b80..000000000 --- a/tests/CurrentDeposition/inputs +++ /dev/null @@ -1,6 +0,0 @@ - -interpolation.nox = 1 -interpolation.noy = 1 -interpolation.noz = 1 - -algo.current_deposition = 3 diff --git a/tests/CurrentDeposition/main.cpp b/tests/CurrentDeposition/main.cpp deleted file mode 100644 index 2b80ce393..000000000 --- a/tests/CurrentDeposition/main.cpp +++ /dev/null @@ -1,139 +0,0 @@ - -#include <random> - -#include <AMReX.H> -#include <AMReX_ParmParse.H> -#include <AMReX_Vector.H> -#include <AMReX_MultiFab.H> -#include <AMReX_MultiFabUtil.H> -#include <AMReX_PlotFileUtil.H> - -#include <WarpX_f.H> -#include <WarpXConst.H> - -using namespace amrex; - -int main(int argc, char* argv[]) -{ - amrex::Initialize(argc,argv); - - { - long nox=1, noy=1, noz=1; - { - ParmParse pp("interpolation"); - pp.query("nox", nox); - pp.query("noy", noy); - pp.query("noz", noz); - if (nox != noy || nox != noz) { - amrex::Abort("warpx.nox, noy and noz must be equal"); - } - if (nox < 1) { - amrex::Abort("warpx.nox must >= 1"); - } - } - - Real charge = -PhysConst::q_e; - Real weight = 10.0; - Real dt = 1.e-10; - - long current_deposition_algo = 3; - { - ParmParse pp("algo"); - pp.query("current_deposition", current_deposition_algo); - } - - long nx = 64, ny = 64, nz = 64; - long np = nx*ny*nz; - - Real xyzmin[3] = {0.5, 1.4, 0.3}; - - Vector<Real> xp, yp, zp, uxp, uyp, uzp, giv, wp; - Real dx[3] = {1.0/nx, 1.0/ny, 1.0/nz}; - - std::mt19937 rand_eng(42); - std::uniform_real_distribution<Real> rand_dis(0.0,1.0); - - for (int k = 0; k < nz; ++k) { - for (int j = 0; j < ny; ++j) { - for (int i = 0; i < nx; ++i) { - Real x0 = xyzmin[0] + i*dx[0]; - Real y0 = xyzmin[1] + j*dx[1]; - Real z0 = xyzmin[2] + k*dx[2]; - xp.push_back(x0 + dx[0]*rand_dis(rand_eng)); - yp.push_back(y0 + dx[1]*rand_dis(rand_eng)); - zp.push_back(z0 + dx[2]*rand_dis(rand_eng)); - wp.push_back(weight); - Real vx,vy,vz,v2; - do { - vx = rand_dis(rand_eng); - vy = rand_dis(rand_eng); - vz = rand_dis(rand_eng); - v2 = vx*vx + vy*vy + vz*vz; - } while(v2 >= 0.999999); - Real gam = 1.0/sqrt(1.0-v2); - uxp.push_back(vx*gam); - uyp.push_back(vy*gam); - uzp.push_back(vz*gam); - giv.push_back(1.0/gam); - } - } - } - - const int ng = nox; - Box domain_box {IntVect{D_DECL(0,0,0)}, IntVect{D_DECL(nx-1,ny-1,nz-1)}}; - -#if (BL_SPACEDIM == 3) - IntVect jx_nodal_flag(0,1,1); - IntVect jy_nodal_flag(1,0,1); - IntVect jz_nodal_flag(1,1,0); -#elif (BL_SPACEDIM == 2) - IntVect jx_nodal_flag(0,1); // x is the first dimension to AMReX - IntVect jy_nodal_flag(1,1); // y is the missing dimension to 2D AMReX - IntVect jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX -#endif - - BoxArray ba{domain_box}; - DistributionMapping dm{ba}; - - MultiFab jx(amrex::convert(ba,jx_nodal_flag), dm, 1, ng); - MultiFab jy(amrex::convert(ba,jy_nodal_flag), dm, 1, ng); - MultiFab jz(amrex::convert(ba,jz_nodal_flag), dm, 1, ng); - - FArrayBox& jxfab = jx[0]; - FArrayBox& jyfab = jy[0]; - FArrayBox& jzfab = jz[0]; - jxfab.setVal(0.0); - jyfab.setVal(0.0); - jzfab.setVal(0.0); - - long ngx = ng; - long ngy = ng; - long ngz = ng; - long lvect = 8; - auto jxntot = jxfab.length(); - auto jyntot = jyfab.length(); - auto jzntot = jzfab.length(); - warpx_current_deposition(jxfab.dataPtr(), &ngx, jxntot.getVect(), - jyfab.dataPtr(), &ngy, jyntot.getVect(), - jzfab.dataPtr(), &ngz, jzntot.getVect(), - &np, xp.data(), yp.data(), zp.data(), - uxp.data(), uyp.data(), uzp.data(), - giv.data(), wp.data(), &charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dt, &dx[0], &dx[1], &dx[2], - &nox, &noy,&noz, - &lvect, ¤t_deposition_algo); - - MultiFab plotmf(ba, dm, 3, 0); - amrex::average_edge_to_cellcenter(plotmf, 0, {&jx, &jy, &jz}); - - RealBox realbox{domain_box, dx, xyzmin}; - int is_per[3] = {0,0,0}; - Geometry geom{domain_box, &realbox, 0, is_per}; - std::string plotname{"plotfiles/plt00000"}; - Vector<std::string> varnames{"jx", "jy", "jz"}; - amrex::WriteSingleLevelPlotfile(plotname, plotmf, varnames, geom, 0.0, 0); - } - - amrex::Finalize(); -} diff --git a/tests/FieldGather/GNUmakefile b/tests/FieldGather/GNUmakefile deleted file mode 100644 index 05ac27a4b..000000000 --- a/tests/FieldGather/GNUmakefile +++ /dev/null @@ -1,25 +0,0 @@ -AMREX_HOME ?= ../../../amrex -PICSAR_HOME ?= ../../../picsar - -USE_PARTICLES = TRUE - -DEBUG = FALSE -USE_MPI = TRUE -USE_OMP = TRUE -PROFILE = FALSE -COMP = gnu -DIM = 3 -PRECISION = DOUBLE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package -include $(AMREX_HOME)/Src/Base/Make.package -include $(PICSAR_HOME)/src/Make.package - -DEFINES += -DWARPX - -default: $(executable) - @echo SUCCESS - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/tests/FieldGather/Make.package b/tests/FieldGather/Make.package deleted file mode 100644 index 5d4d90336..000000000 --- a/tests/FieldGather/Make.package +++ /dev/null @@ -1,9 +0,0 @@ - -CEXE_sources += main.cpp - -CEXE_headers += WarpX_f.H - -F90EXE_sources += WarpX_picsar.F90 - -INCLUDE_LOCATIONS += ../../Source -VPATH_LOCATIONS += ../../Source diff --git a/tests/FieldGather/inputs b/tests/FieldGather/inputs deleted file mode 100644 index 8bd6772ab..000000000 --- a/tests/FieldGather/inputs +++ /dev/null @@ -1,6 +0,0 @@ - -interpolation.nox = 1 -interpolation.noy = 1 -interpolation.noz = 1 - -algo.field_gathering = 1 diff --git a/tests/FieldGather/main.cpp b/tests/FieldGather/main.cpp deleted file mode 100644 index e0719ae9f..000000000 --- a/tests/FieldGather/main.cpp +++ /dev/null @@ -1,198 +0,0 @@ - -#include <random> - -#include <AMReX.H> -#include <AMReX_ParmParse.H> -#include <AMReX_Vector.H> -#include <AMReX_MultiFab.H> -#include <AMReX_BoxIterator.H> -#include <AMReX_PlotFileUtil.H> - -#include <WarpX_f.H> - -using namespace amrex; - -int main(int argc, char* argv[]) -{ - amrex::Initialize(argc,argv); - - { - long nox=1, noy=1, noz=1; - { - ParmParse pp("interpolation"); - pp.query("nox", nox); - pp.query("noy", noy); - pp.query("noz", noz); - if (nox != noy || nox != noz) { - amrex::Abort("warpx.nox, noy and noz must be equal"); - } - if (nox < 1) { - amrex::Abort("warpx.nox must >= 1"); - } - } - - long field_gathering_algo = 1; - - { - ParmParse pp("algo"); - pp.query("field_gathering", field_gathering_algo); - } - - long nx = 64, ny = 64, nz = 64; - long np = nx*ny*nz; - - Real xyzmin[3] = {0.5, 1.4, 0.3}; - int ixyzmin[3] = {0, 0, 0}; - - Vector<Real> xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp; - Real dx[3] = {1.0/nx, 1.0/ny, 1.0/nz}; - - std::mt19937 rand_eng(42); - std::uniform_real_distribution<Real> rand_dis(0.0,dx[0]); - - for (int k = 0; k < nz; ++k) { - for (int j = 0; j < ny; ++j) { - for (int i = 0; i < nx; ++i) { - Real x0 = xyzmin[0] + i*dx[0]; - Real y0 = xyzmin[1] + j*dx[1]; - Real z0 = xyzmin[2] + k*dx[2]; - xp.push_back(x0 + rand_dis(rand_eng)); - yp.push_back(y0 + rand_dis(rand_eng)); - zp.push_back(z0 + rand_dis(rand_eng)); - } - } - } - - Exp.resize(np,0.0); - Eyp.resize(np,0.0); - Ezp.resize(np,0.0); - Bxp.resize(np,0.0); - Byp.resize(np,0.0); - Bzp.resize(np,0.0); - - const int ng = nox; - Box domain_box {IntVect{D_DECL(0,0,0)}, IntVect{D_DECL(nx-1,ny-1,nz-1)}}; - -#if (BL_SPACEDIM == 3) - IntVect Bx_nodal_flag(1,0,0); - IntVect By_nodal_flag(0,1,0); - IntVect Bz_nodal_flag(0,0,1); -#elif (BL_SPACEDIM == 2) - IntVect Bx_nodal_flag(1,0); // x is the first dimension to AMReX - IntVect By_nodal_flag(0,0); // y is the missing dimension to 2D AMReX - IntVect Bz_nodal_flag(0,1); // z is the second dimension to 2D AMReX -#endif - -#if (BL_SPACEDIM == 3) - IntVect Ex_nodal_flag(0,1,1); - IntVect Ey_nodal_flag(1,0,1); - IntVect Ez_nodal_flag(1,1,0); -#elif (BL_SPACEDIM == 2) - IntVect Ex_nodal_flag(0,1); // x is the first dimension to AMReX - IntVect Ey_nodal_flag(1,1); // y is the missing dimension to 2D AMReX - IntVect Ez_nodal_flag(1,0); // z is the second dimension to 2D AMReX -#endif - - BoxArray ba{domain_box}; - DistributionMapping dm{ba}; - - MultiFab Bx(amrex::convert(ba,Bx_nodal_flag), dm, 1, ng); - MultiFab By(amrex::convert(ba,By_nodal_flag), dm, 1, ng); - MultiFab Bz(amrex::convert(ba,Bz_nodal_flag), dm, 1, ng); - MultiFab Ex(amrex::convert(ba,Ex_nodal_flag), dm, 1, ng); - MultiFab Ey(amrex::convert(ba,Ey_nodal_flag), dm, 1, ng); - MultiFab Ez(amrex::convert(ba,Ez_nodal_flag), dm, 1, ng); - - FArrayBox& exfab = Ex[0]; - FArrayBox& eyfab = Ey[0]; - FArrayBox& ezfab = Ez[0]; - FArrayBox& bxfab = Bx[0]; - FArrayBox& byfab = By[0]; - FArrayBox& bzfab = Bz[0]; - - std::uniform_real_distribution<Real> rand_dis2(0.0,100.0); - - for (BoxIterator bxi(exfab.box()); bxi.ok(); ++bxi) - { - exfab(bxi()) = rand_dis2(rand_eng); - } - - for (BoxIterator bxi(eyfab.box()); bxi.ok(); ++bxi) - { - eyfab(bxi()) = rand_dis2(rand_eng); - } - - for (BoxIterator bxi(ezfab.box()); bxi.ok(); ++bxi) - { - ezfab(bxi()) = rand_dis2(rand_eng); - } - - for (BoxIterator bxi(bxfab.box()); bxi.ok(); ++bxi) - { - bxfab(bxi()) = rand_dis2(rand_eng); - } - - for (BoxIterator bxi(byfab.box()); bxi.ok(); ++bxi) - { - byfab(bxi()) = rand_dis2(rand_eng); - } - - for (BoxIterator bxi(bzfab.box()); bxi.ok(); ++bxi) - { - bzfab(bxi()) = rand_dis2(rand_eng); - } - - long ngx = ng; - long ngy = ng; - long ngz = ng; - - const int ll4symtry = false; - const int l_lower_order_in_v = true; - const int l_nodal = false; - long lvect_fieldgathe = 64; - warpx_geteb_energy_conserving(&np, xp.data(), yp.data(), zp.data(), - Exp.data(),Eyp.data(),Ezp.data(), - Bxp.data(),Byp.data(),Bzp.data(), - ixyzmin, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], - &nox, &noy, &noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &l_lower_order_in_v, - &l_nodal, - &lvect_fieldgathe, - &field_gathering_algo); - - MultiFab plotmf(ba, dm, 6, 0); - FArrayBox& plotfab = plotmf[0]; - plotfab.setVal(0.0); - for (int k = 0; k < nz; ++k) { - for (int j = 0; j < ny; ++j) { - for (int i = 0; i < nx; ++i) { - IntVect cell{i,j,k}; - int ip = i + j*nx + k*nx*ny; - plotfab(cell,0) += Exp[ip]; - plotfab(cell,1) += Eyp[ip]; - plotfab(cell,2) += Ezp[ip]; - plotfab(cell,3) += Bxp[ip]; - plotfab(cell,4) += Byp[ip]; - plotfab(cell,5) += Bzp[ip]; - } - } - } - - RealBox realbox{domain_box, dx, xyzmin}; - int is_per[3] = {0,0,0}; - Geometry geom{domain_box, &realbox, 0, is_per}; - std::string plotname{"plotfiles/plt00000"}; - Vector<std::string> varnames{"Ex", "Ey", "Ez", "Bx", "By", "Bz"}; - amrex::WriteSingleLevelPlotfile(plotname, plotmf, varnames, geom, 0.0, 0); - } - - amrex::Finalize(); -} diff --git a/tests/FieldSolver/GNUmakefile b/tests/FieldSolver/GNUmakefile deleted file mode 100644 index 05ac27a4b..000000000 --- a/tests/FieldSolver/GNUmakefile +++ /dev/null @@ -1,25 +0,0 @@ -AMREX_HOME ?= ../../../amrex -PICSAR_HOME ?= ../../../picsar - -USE_PARTICLES = TRUE - -DEBUG = FALSE -USE_MPI = TRUE -USE_OMP = TRUE -PROFILE = FALSE -COMP = gnu -DIM = 3 -PRECISION = DOUBLE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package -include $(AMREX_HOME)/Src/Base/Make.package -include $(PICSAR_HOME)/src/Make.package - -DEFINES += -DWARPX - -default: $(executable) - @echo SUCCESS - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/tests/FieldSolver/Make.package b/tests/FieldSolver/Make.package deleted file mode 100644 index 3129dcc37..000000000 --- a/tests/FieldSolver/Make.package +++ /dev/null @@ -1,9 +0,0 @@ - -CEXE_sources += main.cpp - -CEXE_headers += WarpX_f.H WarpXConst.H - -F90EXE_sources += WarpX_picsar.F90 - -INCLUDE_LOCATIONS += ../../Source -VPATH_LOCATIONS += ../../Source diff --git a/tests/FieldSolver/inputs b/tests/FieldSolver/inputs deleted file mode 100644 index 7e8cf4580..000000000 --- a/tests/FieldSolver/inputs +++ /dev/null @@ -1,4 +0,0 @@ - -interpolation.nox = 1 -interpolation.noy = 1 -interpolation.noz = 1 diff --git a/tests/FieldSolver/main.cpp b/tests/FieldSolver/main.cpp deleted file mode 100644 index c547d2631..000000000 --- a/tests/FieldSolver/main.cpp +++ /dev/null @@ -1,229 +0,0 @@ - -#include <limits> -#include <random> - -#include <AMReX.H> -#include <AMReX_ParmParse.H> -#include <AMReX_Vector.H> -#include <AMReX_MultiFab.H> -#include <AMReX_MultiFabUtil.H> -#include <AMReX_PlotFileUtil.H> - -#include <WarpXConst.H> -#include <WarpX_f.H> - -using namespace amrex; - -int main(int argc, char* argv[]) -{ - amrex::Initialize(argc,argv); - - { - long nox=1, noy=1, noz=1; - { - ParmParse pp("interpolation"); - pp.query("nox", nox); - pp.query("noy", noy); - pp.query("noz", noz); - if (nox != noy || nox != noz) { - amrex::Abort("warpx.nox, noy and noz must be equal"); - } - if (nox < 1) { - amrex::Abort("warpx.nox must >= 1"); - } - } - - std::mt19937 rand_eng(42); - std::uniform_real_distribution<Real> rand_dis(0.0,1.0); - - long nx = 64, ny = 64, nz = 64; - Real dx[3] = {1.0/nx, 1.0/ny, 1.0/nz}; - Real dt = 1.e-6; - - Box cc_domain{IntVect{D_DECL(0,0,0)},IntVect{D_DECL(nx-1,ny-1,nz-1)}}; - BoxArray grids{cc_domain}; - grids.maxSize(32); - - DistributionMapping dmap {grids}; - - const int ng = nox; - -#if (BL_SPACEDIM == 3) - IntVect Bx_nodal_flag(1,0,0); - IntVect By_nodal_flag(0,1,0); - IntVect Bz_nodal_flag(0,0,1); -#elif (BL_SPACEDIM == 2) - IntVect Bx_nodal_flag(1,0); // x is the first dimension to AMReX - IntVect By_nodal_flag(0,0); // y is the missing dimension to 2D AMReX - IntVect Bz_nodal_flag(0,1); // z is the second dimension to 2D AMReX -#endif - -#if (BL_SPACEDIM == 3) - IntVect Ex_nodal_flag(0,1,1); - IntVect Ey_nodal_flag(1,0,1); - IntVect Ez_nodal_flag(1,1,0); -#elif (BL_SPACEDIM == 2) - IntVect Ex_nodal_flag(0,1); // x is the first dimension to AMReX - IntVect Ey_nodal_flag(1,1); // y is the missing dimension to 2D AMReX - IntVect Ez_nodal_flag(1,0); // z is the second dimension to 2D AMReX -#endif - -#if (BL_SPACEDIM == 3) - IntVect jx_nodal_flag(0,1,1); - IntVect jy_nodal_flag(1,0,1); - IntVect jz_nodal_flag(1,1,0); -#elif (BL_SPACEDIM == 2) - IntVect jx_nodal_flag(0,1); // x is the first dimension to AMReX - IntVect jy_nodal_flag(1,1); // y is the missing dimension to 2D AMReX - IntVect jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX -#endif - - Vector<std::unique_ptr<MultiFab> > current(3); - Vector<std::unique_ptr<MultiFab> > Efield(3); - Vector<std::unique_ptr<MultiFab> > Bfield(3); - // Create the MultiFabs for B - Bfield[0].reset( new MultiFab(amrex::convert(grids,Bx_nodal_flag),dmap,1,ng)); - Bfield[1].reset( new MultiFab(amrex::convert(grids,By_nodal_flag),dmap,1,ng)); - Bfield[2].reset( new MultiFab(amrex::convert(grids,Bz_nodal_flag),dmap,1,ng)); - - // Create the MultiFabs for E - Efield[0].reset( new MultiFab(amrex::convert(grids,Ex_nodal_flag),dmap,1,ng)); - Efield[1].reset( new MultiFab(amrex::convert(grids,Ey_nodal_flag),dmap,1,ng)); - Efield[2].reset( new MultiFab(amrex::convert(grids,Ez_nodal_flag),dmap,1,ng)); - - // Create the MultiFabs for the current - current[0].reset( new MultiFab(amrex::convert(grids,jx_nodal_flag),dmap,1,ng)); - current[1].reset( new MultiFab(amrex::convert(grids,jy_nodal_flag),dmap,1,ng)); - current[2].reset( new MultiFab(amrex::convert(grids,jz_nodal_flag),dmap,1,ng)); - - Bfield[0]->setVal(0.0); - Bfield[1]->setVal(0.0); - Bfield[2]->setVal(0.0); - Efield[0]->setVal(0.0); - Efield[1]->setVal(0.0); - Efield[2]->setVal(0.0); - current[0]->setVal(0.0); - current[1]->setVal(0.0); - current[2]->setVal(0.0); - - for (MFIter mfi(*current[0]); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - FArrayBox& exfab = (*Efield [0])[mfi]; - FArrayBox& eyfab = (*Efield [1])[mfi]; - FArrayBox& ezfab = (*Efield [2])[mfi]; - FArrayBox& bxfab = (*Bfield [0])[mfi]; - FArrayBox& byfab = (*Bfield [1])[mfi]; - FArrayBox& bzfab = (*Bfield [2])[mfi]; - FArrayBox& jxfab = (*current[0])[mfi]; - FArrayBox& jyfab = (*current[1])[mfi]; - FArrayBox& jzfab = (*current[2])[mfi]; - for (IntVect cell=bx.smallEnd(); cell <= bx.bigEnd(); bx.next(cell)) - { - exfab(cell) = rand_dis(rand_eng)*1.e5; - eyfab(cell) = rand_dis(rand_eng)*1.e5; - ezfab(cell) = rand_dis(rand_eng)*1.e5; - bxfab(cell) = rand_dis(rand_eng)*1.e-5; - byfab(cell) = rand_dis(rand_eng)*1.e-5; - bzfab(cell) = rand_dis(rand_eng)*1.e-5; - jxfab(cell) = rand_dis(rand_eng)*1.e10; - jyfab(cell) = rand_dis(rand_eng)*1.e10; - jzfab(cell) = rand_dis(rand_eng)*1.e10; - } - } - - - { // Evolve B - Real dtsdx[3]; -#if (BL_SPACEDIM == 3) - dtsdx[0] = dt / dx[0]; - dtsdx[1] = dt / dx[1]; - dtsdx[2] = dt / dx[2]; -#elif (BL_SPACEDIM == 2) - dtsdx[0] = dt / dx[0]; - dtsdx[1] = std::numeric_limits<Real>::quiet_NaN(); - dtsdx[2] = dt / dx[1]; -#endif - - const int norder = 2; - // FDT Yee solver - const int maxwell_fdtd_solver_id = 0; - - for ( MFIter mfi(*Bfield[0],true); mfi.isValid(); ++mfi ) - { - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - // Call picsar routine for each tile - warpx_push_bvec( - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*Efield[0])[mfi]), - BL_TO_FORTRAN_3D((*Efield[1])[mfi]), - BL_TO_FORTRAN_3D((*Efield[2])[mfi]), - BL_TO_FORTRAN_3D((*Bfield[0])[mfi]), - BL_TO_FORTRAN_3D((*Bfield[1])[mfi]), - BL_TO_FORTRAN_3D((*Bfield[2])[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &maxwell_fdtd_solver_id); - } - } - - { // evolve E - Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; - - Real dtsdx_c2[3]; -#if (BL_SPACEDIM == 3) - dtsdx_c2[0] = (PhysConst::c*PhysConst::c) * dt / dx[0]; - dtsdx_c2[1] = (PhysConst::c*PhysConst::c) * dt / dx[1]; - dtsdx_c2[2] = (PhysConst::c*PhysConst::c) * dt / dx[2]; -#else - dtsdx_c2[0] = (PhysConst::c*PhysConst::c) * dt / dx[0]; - dtsdx_c2[1] = std::numeric_limits<Real>::quiet_NaN(); - dtsdx_c2[2] = (PhysConst::c*PhysConst::c) * dt / dx[1]; -#endif - - const int norder = 2; - - for ( MFIter mfi(*Efield[0],true); mfi.isValid(); ++mfi ) - { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - - // Call picsar routine for each tile - warpx_push_evec( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Efield[0])[mfi]), - BL_TO_FORTRAN_3D((*Efield[1])[mfi]), - BL_TO_FORTRAN_3D((*Efield[2])[mfi]), - BL_TO_FORTRAN_3D((*Bfield[0])[mfi]), - BL_TO_FORTRAN_3D((*Bfield[1])[mfi]), - BL_TO_FORTRAN_3D((*Bfield[2])[mfi]), - BL_TO_FORTRAN_3D((*current[0])[mfi]), - BL_TO_FORTRAN_3D((*current[1])[mfi]), - BL_TO_FORTRAN_3D((*current[2])[mfi]), - &mu_c2_dt, - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - } - } - - MultiFab plotmf(grids, dmap, 6, 0); - amrex::average_edge_to_cellcenter(plotmf, 0, amrex::GetVecOfConstPtrs(Efield)); - amrex::average_face_to_cellcenter(plotmf, 3, amrex::GetVecOfConstPtrs(Bfield)); - - Real xyzmin[3] = {0.0,0.0,0.0}; - RealBox realbox{cc_domain, dx, xyzmin}; - int is_per[3] = {0,0,0}; - Geometry geom{cc_domain, &realbox, 0, is_per}; - std::string plotname{"plotfiles/plt00000"}; - Vector<std::string> varnames{"Ex", "Ey", "Ez", "Bx", "By", "Bz"}; - amrex::WriteSingleLevelPlotfile(plotname, plotmf, varnames, geom, 0.0, 0); - } - - amrex::Finalize(); -} diff --git a/tests/ParticlePusher/GNUmakefile b/tests/ParticlePusher/GNUmakefile deleted file mode 100644 index 05ac27a4b..000000000 --- a/tests/ParticlePusher/GNUmakefile +++ /dev/null @@ -1,25 +0,0 @@ -AMREX_HOME ?= ../../../amrex -PICSAR_HOME ?= ../../../picsar - -USE_PARTICLES = TRUE - -DEBUG = FALSE -USE_MPI = TRUE -USE_OMP = TRUE -PROFILE = FALSE -COMP = gnu -DIM = 3 -PRECISION = DOUBLE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package -include $(AMREX_HOME)/Src/Base/Make.package -include $(PICSAR_HOME)/src/Make.package - -DEFINES += -DWARPX - -default: $(executable) - @echo SUCCESS - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/tests/ParticlePusher/Make.package b/tests/ParticlePusher/Make.package deleted file mode 100644 index 5d4d90336..000000000 --- a/tests/ParticlePusher/Make.package +++ /dev/null @@ -1,9 +0,0 @@ - -CEXE_sources += main.cpp - -CEXE_headers += WarpX_f.H - -F90EXE_sources += WarpX_picsar.F90 - -INCLUDE_LOCATIONS += ../../Source -VPATH_LOCATIONS += ../../Source diff --git a/tests/ParticlePusher/inputs b/tests/ParticlePusher/inputs deleted file mode 100644 index 675def443..000000000 --- a/tests/ParticlePusher/inputs +++ /dev/null @@ -1,2 +0,0 @@ - -algo.particle_pusher = 0 diff --git a/tests/ParticlePusher/main.cpp b/tests/ParticlePusher/main.cpp deleted file mode 100644 index 4df85b81c..000000000 --- a/tests/ParticlePusher/main.cpp +++ /dev/null @@ -1,119 +0,0 @@ - -#include <limits> -#include <random> - -#include <AMReX.H> -#include <AMReX_ParmParse.H> -#include <AMReX_Vector.H> -#include <AMReX_MultiFab.H> -#include <AMReX_PlotFileUtil.H> - -#include <WarpXConst.H> -#include <WarpX_f.H> - -using namespace amrex; - -int main(int argc, char* argv[]) -{ - amrex::Initialize(argc,argv); - - { - long particle_pusher_algo = 0; - { - ParmParse pp("algo"); - pp.query("particle_pusher", particle_pusher_algo); - } - - long nx = 64, ny = 64, nz = 64; - long np = nx*ny*nz; - - Real xyzmin[3] = {0.5, 1.4, 0.3}; - - Vector<Real> xp, yp, zp, uxp, uyp, uzp, giv, Exp, Eyp, Ezp, Bxp, Byp, Bzp; - Real dx[3] = {1.0/nx, 1.0/ny, 1.0/nz}; - - std::mt19937 rand_eng(42); - std::uniform_real_distribution<Real> rand_dis(0.0,1.0); - - for (int k = 0; k < nz; ++k) { - for (int j = 0; j < ny; ++j) { - for (int i = 0; i < nx; ++i) { - Real x0 = xyzmin[0] + i*dx[0]; - Real y0 = xyzmin[1] + j*dx[1]; - Real z0 = xyzmin[2] + k*dx[2]; - xp.push_back(x0 + dx[0]*rand_dis(rand_eng)); - yp.push_back(y0 + dx[1]*rand_dis(rand_eng)); - zp.push_back(z0 + dx[2]*rand_dis(rand_eng)); - Real vx,vy,vz,v2; - do { - vx = rand_dis(rand_eng); - vy = rand_dis(rand_eng); - vz = rand_dis(rand_eng); - v2 = vx*vx + vy*vy + vz*vz; - } while(v2 >= 0.999999); - Real gam = 1.0/sqrt(1.0-v2); - uxp.push_back(vx*PhysConst::c*gam); - uyp.push_back(vy*PhysConst::c*gam); - uzp.push_back(vz*PhysConst::c*gam); - giv.push_back(std::numeric_limits<Real>::quiet_NaN()); - } - } - } - - Exp.resize(np,0.0); - Eyp.resize(np,0.0); - Ezp.resize(np,0.0); - Bxp.resize(np,0.0); - Byp.resize(np,0.0); - Bzp.resize(np,0.0); - for (int i = 0; i < np; ++i) { - Exp[i] = rand_dis(rand_eng)*1.0e5; - Eyp[i] = rand_dis(rand_eng)*1.0e5; - Ezp[i] = rand_dis(rand_eng)*1.0e5; - Bxp[i] = rand_dis(rand_eng)*1.0e-5; - Byp[i] = rand_dis(rand_eng)*1.0e-5; - Bzp[i] = rand_dis(rand_eng)*1.0e-5; - } - - Real charge = -PhysConst::q_e; - Real mass = PhysConst::m_e; - Real dt = 1.e-10; - warpx_particle_pusher(&np, xp.data(), yp.data(), zp.data(), - uxp.data(), uyp.data(), uzp.data(), giv.data(), - Exp.dataPtr(), Eyp.dataPtr(), Ezp.dataPtr(), - Bxp.dataPtr(), Byp.dataPtr(), Bzp.dataPtr(), - &charge, &mass, &dt, - &particle_pusher_algo); - - Box plotbox{IntVect{D_DECL(0,0,0)},IntVect{D_DECL(nx-1,ny-1,nz-1)}}; - BoxArray plotba {plotbox}; - DistributionMapping plotdm {plotba}; - MultiFab plotmf(plotba, plotdm, 7, 0); - FArrayBox& plotfab = plotmf[0]; - plotfab.setVal(0.0); - for (int k = 0; k < nz; ++k) { - for (int j = 0; j < ny; ++j) { - for (int i = 0; i < nx; ++i) { - IntVect cell{i,j,k}; - int ip = i + j*nx + k*nx*ny; - plotfab(cell,0) += xp[ip]; - plotfab(cell,1) += yp[ip]; - plotfab(cell,2) += zp[ip]; - plotfab(cell,3) += uxp[ip]; - plotfab(cell,4) += uyp[ip]; - plotfab(cell,5) += uzp[ip]; - plotfab(cell,6) += 1.0/giv[ip]; - } - } - } - - RealBox realbox{plotbox, dx, xyzmin}; - int is_per[3] = {0,0,0}; - Geometry geom{plotbox, &realbox, 0, is_per}; - std::string plotname{"plotfiles/plt00000"}; - Vector<std::string> varnames{"x", "y", "z", "ux", "uy", "uz", "gamma"}; - amrex::WriteSingleLevelPlotfile(plotname, plotmf, varnames, geom, 0.0, 0); - } - - amrex::Finalize(); -} |