From e133d2202685d6e478b42f9b554b00d5fa722801 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 25 Oct 2019 11:56:10 -0700 Subject: replace 'boosted frame diags' with 'back-transformed diags' --- Python/pywarpx/picmi.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'Python/pywarpx/picmi.py') diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index c33700278..c4e6803d5 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -671,10 +671,10 @@ class LabFrameFieldDiagnostic(picmistandard.PICMI_LabFrameFieldDiagnostic): pywarpx.warpx.check_consistency('dt_snapshots_lab', self.dt_snapshots, 'The time between snapshots must be the same in all lab frame diagnostics') pywarpx.warpx.check_consistency('lab_data_directory', self.write_dir, 'The write directory must be the same in all lab frame diagnostics') - pywarpx.warpx.do_boosted_frame_diagnostic = 1 + pywarpx.warpx.do_back_transformed_diagnostics = 1 pywarpx.warpx.num_snapshots_lab = self.num_snapshots pywarpx.warpx.dt_snapshots_lab = self.dt_snapshots - pywarpx.warpx.do_boosted_frame_fields = 1 + pywarpx.warpx.do_back_transformed_fields = 1 pywarpx.warpx.lab_data_directory = self.write_dir @@ -685,8 +685,8 @@ class LabFrameParticleDiagnostic(picmistandard.PICMI_LabFrameParticleDiagnostic) pywarpx.warpx.check_consistency('dt_snapshots_lab', self.dt_snapshots, 'The time between snapshots must be the same in all lab frame diagnostics') pywarpx.warpx.check_consistency('lab_data_directory', self.write_dir, 'The write directory must be the same in all lab frame diagnostics') - pywarpx.warpx.do_boosted_frame_diagnostic = 1 + pywarpx.warpx.do_back_transformed_diagnostics = 1 pywarpx.warpx.num_snapshots_lab = self.num_snapshots pywarpx.warpx.dt_snapshots_lab = self.dt_snapshots - pywarpx.warpx.do_boosted_frame_particles = 1 + pywarpx.warpx.do_back_transformed_particles = 1 pywarpx.warpx.lab_data_directory = self.write_dir -- cgit v1.2.3 From bdef308e6c5b4aeed8190b6ecdb25b00a51ca5f9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 6 Nov 2019 16:22:41 -0800 Subject: Implement space-charge initialization (non-relativistic, single-level) --- Docs/source/running_cpp/parameters.rst | 4 + .../space_charge_initialization/analysis.py | 85 +++++++++++++++ .../Modules/space_charge_initialization/inputs | 29 ++++++ Python/pywarpx/picmi.py | 12 +-- Regression/WarpX-tests.ini | 34 ++++++ Source/Initialization/InitSpaceChargeField.cpp | 116 +++++++++++++++++++++ Source/Initialization/Make.package | 2 + Source/Initialization/WarpXInitData.cpp | 83 ++------------- Source/Make.WarpX | 1 + Source/Particles/PhysicalParticleContainer.cpp | 1 + Source/Particles/WarpXParticleContainer.H | 3 +- Source/Particles/WarpXParticleContainer.cpp | 82 ++++++--------- Source/WarpX.H | 22 +--- 13 files changed, 326 insertions(+), 148 deletions(-) create mode 100755 Examples/Modules/space_charge_initialization/analysis.py create mode 100644 Examples/Modules/space_charge_initialization/inputs create mode 100644 Source/Initialization/InitSpaceChargeField.cpp (limited to 'Python/pywarpx/picmi.py') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index e5727e376..50803560d 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -224,6 +224,10 @@ Particle initialization initialization. This can be required whith a moving window and/or when running in a boosted frame. +* ``.initialize_self_fields`` (`0` or `1`) + Whether to calculate the space-charge fields associated with this species + at the beginning of the simulation. + * ``.profile`` (`string`) Density profile for this species. The options are: diff --git a/Examples/Modules/space_charge_initialization/analysis.py b/Examples/Modules/space_charge_initialization/analysis.py new file mode 100755 index 000000000..30954c41d --- /dev/null +++ b/Examples/Modules/space_charge_initialization/analysis.py @@ -0,0 +1,85 @@ +#! /usr/bin/env python +""" +This script checks the space-charge initialization routine, by +verifying that the space-charge field of a Gaussian beam corresponds to +the expected theoretical field. +""" +import sys +import yt +import numpy as np +import matplotlib.pyplot as plt +import scipy.constants as scc +from scipy.special import gammainc +yt.funcs.mylog.setLevel(0) + +# Parameters from the Simulation +Qtot = -1.e-20 +r0 = 2.e-6 + +# Open data file +filename = sys.argv[1] +ds = yt.load( filename ) +# Extract data +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray().squeeze() +if ds.dimensionality == 2: + # Rename the z dimension as y, so as to make this script work for 2d and 3d + Ey_array = ad0['Ez'].to_ndarray().squeeze() +elif ds.dimensionality == 3: + Ey_array = ad0['Ey'].to_ndarray() + Ez_array = ad0['Ez'].to_ndarray() + +# Extract grid coordinates +Nx, Ny, Nz = ds.domain_dimensions +xmin, ymin, zmin = ds.domain_left_edge.v +Lx, Ly, Lz = ds.domain_width.v +x = xmin + Lx/Nx*(0.5+np.arange(Nx)) +y = ymin + Ly/Ny*(0.5+np.arange(Ny)) +z = zmin + Lz/Nz*(0.5+np.arange(Nz)) + +# Compute theoretical field +if ds.dimensionality == 2: + x_2d, y_2d = np.meshgrid(x, y, indexing='ij') + r2 = x_2d**2 + y_2d**2 + factor = (Qtot/r0)/(2*np.pi*scc.epsilon_0*r2) * (1-np.exp(-r2/(2*r0**2))) + Ex_th = x_2d * factor + Ey_th = y_2d * factor +elif ds.dimensionality == 3: + x_2d, y_2d, z_2d = np.meshgrid(x, y, z, indexing='ij') + r2 = x_2d**2 + y_2d**2 + z_2d**2 + factor = Qtot/(4*np.pi*scc.epsilon_0*r2**1.5) * gammainc(3./2, r2/(2.*r0**2)) + Ex_th = factor*x_2d + Ey_th = factor*y_2d + Ez_th = factor*z_2d + +# Plot theory and data +def make_2d(arr): + if arr.ndim == 3: + return arr[:,:,Nz//2] + else: + return arr +plt.figure(figsize=(10,10)) +plt.subplot(221) +plt.title('Ex: Theory') +plt.imshow(make_2d(Ex_th)) +plt.colorbar() +plt.subplot(222) +plt.title('Ex: Simulation') +plt.imshow(make_2d(Ex_array)) +plt.colorbar() +plt.subplot(223) +plt.title('Ey: Theory') +plt.imshow(make_2d(Ey_th)) +plt.colorbar() +plt.subplot(224) +plt.title('Ey: Simulation') +plt.imshow(make_2d(Ey_array)) +plt.colorbar() +plt.savefig('Comparison.png') + +# Automatically check the results +print( 'Relative error in Ex: %.3f'%(abs(Ex_array-Ex_th).max()/Ex_th.max())) +assert np.allclose( Ex_array, Ex_th, atol=0.15*Ex_th.max() ) +assert np.allclose( Ey_array, Ey_th, atol=0.15*Ey_th.max() ) +if ds.dimensionality == 3: + assert np.allclose( Ez_array, Ez_th, atol=0.15*Ez_th.max() ) diff --git a/Examples/Modules/space_charge_initialization/inputs b/Examples/Modules/space_charge_initialization/inputs new file mode 100644 index 000000000..f26f0bf22 --- /dev/null +++ b/Examples/Modules/space_charge_initialization/inputs @@ -0,0 +1,29 @@ +max_step = 0 +amr.n_cell = 64 64 64 +amr.max_grid_size = 32 +amr.max_level = 0 + +amr.plot_int = 1 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 0 # Is periodic? +geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6 +geometry.prob_hi = 50.e-6 50.e-6 50.e-6 + +particles.nspecies = 1 +particles.species_names = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.initialize_self_fields = 1 +beam.x_rms = 2.e-6 +beam.y_rms = 2.e-6 +beam.z_rms = 2.e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = 0.e-6 +beam.npart = 2000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 0.0 diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index c4e6803d5..6e9c9153b 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -56,7 +56,7 @@ class Species(picmistandard.PICMI_Species): if self.mass is None: self.mass = element.mass*periodictable.constants.atomic_mass_constant - def initialize_inputs(self, layout): + def initialize_inputs(self, layout, initialize_self_fields=False): self.species_number = pywarpx.particles.nspecies pywarpx.particles.nspecies += 1 @@ -68,7 +68,8 @@ class Species(picmistandard.PICMI_Species): else: pywarpx.particles.species_names += ' ' + self.name - self.species = pywarpx.Bucket.Bucket(self.name, mass=self.mass, charge=self.charge, injection_style = 'python') + self.species = pywarpx.Bucket.Bucket(self.name, mass=self.mass, charge=self.charge, + injection_style = 'python', initialize_self_fields=int(initialize_self_fields)) pywarpx.Particles.particles_list.append(self.species) if self.initial_distribution is not None: @@ -77,9 +78,9 @@ class Species(picmistandard.PICMI_Species): picmistandard.PICMI_MultiSpecies.Species_class = Species class MultiSpecies(picmistandard.PICMI_MultiSpecies): - def initialize_inputs(self, layout): + def initialize_inputs(self, layout, initialize_self_fields=False): for species in self.species_instances_list: - species.initialize_inputs(layout) + species.initialize_inputs(layout, initialize_self_fields) class GaussianBunchDistribution(picmistandard.PICMI_GaussianBunchDistribution): @@ -544,8 +545,7 @@ class Simulation(picmistandard.PICMI_Simulation): self.solver.initialize_inputs() for i in range(len(self.species)): - assert not self.initialize_self_fields[i], Exception('WarpX does not support initializing self fields') - self.species[i].initialize_inputs(self.layouts[i]) + self.species[i].initialize_inputs(self.layouts[i], self.initialize_self_fields[i]) for i in range(len(self.lasers)): self.lasers[i].initialize_inputs() diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 19da4ba04..c23a751a8 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -616,6 +616,40 @@ compareParticles = 0 particleTypes = electrons tolerance = 1.e-14 +[space_charge_initialization_2d] +buildDir = . +inputFile = Examples/Modules/space_charge_initialization/inputs +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +runtime_params = warpx.do_dynamic_scheduling=0 +analysisRoutine = Examples/Modules/space_charge_initialization/analysis.py +analysisOutputImage = Comparison.png + +[space_charge_initialization] +buildDir = . +inputFile = Examples/Modules/space_charge_initialization/inputs +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +runtime_params = warpx.do_dynamic_scheduling=0 +analysisRoutine = Examples/Modules/space_charge_initialization/analysis.py +analysisOutputImage = Comparison.png + [particles_in_pml_2d] buildDir = . inputFile = Examples/Tests/particles_in_PML/inputs.2d diff --git a/Source/Initialization/InitSpaceChargeField.cpp b/Source/Initialization/InitSpaceChargeField.cpp new file mode 100644 index 000000000..771078af3 --- /dev/null +++ b/Source/Initialization/InitSpaceChargeField.cpp @@ -0,0 +1,116 @@ + +#include +#include +#include + +#include + +using namespace amrex; + +void +WarpX::InitSpaceChargeField (WarpXParticleContainer& pc) +{ + // Allocate fields for charge and potential + const int num_levels = max_level + 1; + Vector > rho(num_levels); + Vector > phi(num_levels); + const int ng = WarpX::nox; + for (int lev = 0; lev <= max_level; lev++) { + BoxArray nba = boxArray(lev); + nba.surroundingNodes(); + rho[lev].reset(new MultiFab(nba, dmap[lev], 1, ng)); // Make ng big enough/use rho from sim + phi[lev].reset(new MultiFab(nba, dmap[lev], 1, 0)); + phi[lev]->setVal(0.); + } + + // Deposit particle charge density (source of Poisson solver) + bool local = false; + bool reset = true; + pc.DepositCharge(rho, local, reset); + + // Call amrex's multigrid solver + // ----------------------------- + + // Define the boundary conditions + Array lobc, hibc; + for (int idim=0; idimarray(mfi); + const auto& Ex_arr = (*Efield_fp[lev][0])[mfi].array(); + const auto& Ey_arr = (*Efield_fp[lev][1])[mfi].array(); + const auto& Ez_arr = (*Efield_fp[lev][2])[mfi].array(); + +#if (AMREX_SPACEDIM == 3) + amrex::ParallelFor( tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ex_arr(i,j,k) += inv_dx*( phi_arr(i+1,j,k) - phi_arr(i,j,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ey_arr(i,j,k) += inv_dy*( phi_arr(i,j+1,k) - phi_arr(i,j,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ez_arr(i,j,k) += inv_dz*( phi_arr(i,j,k+1) - phi_arr(i,j,k) ); + } + ); +#else + amrex::ParallelFor( tbx, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ex_arr(i,j,k) += inv_dx*( phi_arr(i+1,j,k) - phi_arr(i,j,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ez_arr(i,j,k) += inv_dz*( phi_arr(i,j+1,k) - phi_arr(i,j,k) ); + } + ); +#endif + } + } +} diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 2c6458b6d..0b65c1ab1 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -14,5 +14,7 @@ CEXE_sources += InjectorMomentum.cpp CEXE_headers += CustomDensityProb.H CEXE_headers += CustomMomentumProb.H +CEXE_sources += InitSpaceChargeField.cpp + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Initialization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Initialization diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 0814f369b..4fd5d066d 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -111,9 +111,13 @@ WarpX::InitFromScratch () mypc->AllocData(); mypc->InitData(); -#ifdef USE_OPENBC_POISSON - InitOpenbc(); -#endif + // Loop through species and calculate their space-charge field + for (int ispecies=0; ispeciesnSpecies(); ispecies++){ + WarpXParticleContainer& species = mypc->GetParticleContainer(ispecies); + if (species.initialize_self_fields) { + InitSpaceChargeField(species); + } + } InitPML(); @@ -226,79 +230,6 @@ WarpX::PostRestart () mypc->PostRestart(); } -#ifdef USE_OPENBC_POISSON -void -WarpX::InitOpenbc () -{ -#ifndef BL_USE_MPI - static_assert(false, "must use MPI"); -#endif - - static_assert(AMREX_SPACEDIM == 3, "Openbc is 3D only"); - BL_ASSERT(finestLevel() == 0); - - const int lev = 0; - - const Geometry& gm = Geom(lev); - const Box& gbox = gm.Domain(); - int lohi[6]; - warpx_openbc_decompose(gbox.loVect(), gbox.hiVect(), lohi, lohi+3); - - int nprocs = ParallelDescriptor::NProcs(); - int myproc = ParallelDescriptor::MyProc(); - Vector alllohi(6*nprocs,100000); - - MPI_Allgather(lohi, 6, MPI_INT, alllohi.data(), 6, MPI_INT, ParallelDescriptor::Communicator()); - - BoxList bl{IndexType::TheNodeType()}; - for (int i = 0; i < nprocs; ++i) - { - bl.push_back(Box(IntVect(alllohi[6*i ],alllohi[6*i+1],alllohi[6*i+2]), - IntVect(alllohi[6*i+3],alllohi[6*i+4],alllohi[6*i+5]), - IndexType::TheNodeType())); - } - BoxArray ba{bl}; - - Vector iprocmap(nprocs+1); - std::iota(iprocmap.begin(), iprocmap.end(), 0); - iprocmap.back() = myproc; - - DistributionMapping dm{iprocmap}; - - MultiFab rho_openbc(ba, dm, 1, 0); - MultiFab phi_openbc(ba, dm, 1, 0); - - bool local = true; - const std::unique_ptr& rho = mypc->GetChargeDensity(lev, local); - - rho_openbc.setVal(0.0); - rho_openbc.copy(*rho, 0, 0, 1, rho->nGrow(), 0, gm.periodicity(), FabArrayBase::ADD); - - const Real* dx = gm.CellSize(); - - warpx_openbc_potential(rho_openbc[myproc].dataPtr(), phi_openbc[myproc].dataPtr(), dx); - - BoxArray nba = boxArray(lev); - nba.surroundingNodes(); - MultiFab phi(nba, DistributionMap(lev), 1, 0); - phi.copy(phi_openbc, gm.periodicity()); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(phi); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - warpx_compute_E(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_3D(phi[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][0])[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][1])[mfi]), - BL_TO_FORTRAN_3D((*Efield[lev][2])[mfi]), - dx); - } -} -#endif - void WarpX::InitLevelData (int lev, Real time) { diff --git a/Source/Make.WarpX b/Source/Make.WarpX index b81fed147..ec004f8fd 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -76,6 +76,7 @@ include $(AMREX_HOME)/Src/Base/Make.package include $(AMREX_HOME)/Src/Particle/Make.package include $(AMREX_HOME)/Src/Boundary/Make.package include $(AMREX_HOME)/Src/AmrCore/Make.package +include $(AMREX_HOME)/Src/LinearSolvers/MLMG/Make.package ifeq ($(USE_SENSEI_INSITU),TRUE) include $(AMREX_HOME)/Src/Amr/Make.package diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 51690d659..6c2aef8a6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -41,6 +41,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); + pp.query("initialize_self_fields", initialize_self_fields); // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index dbd913c5b..cdfe112af 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -185,7 +185,7 @@ public: const amrex::MultiFab& Bz) = 0; void DepositCharge(amrex::Vector >& rho, - bool local = false); + bool local = false, bool reset = true); std::unique_ptr GetChargeDensity(int lev, bool local = false); virtual void DepositCharge(WarpXParIter& pti, @@ -254,6 +254,7 @@ public: void setNextID(int next_id) { ParticleType::NextID(next_id); } bool do_splitting = false; + bool initialize_self_fields = false; // split along diagonals (0) or axes (1) int split_type = 0; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 7fb57500d..29a7d95c7 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -501,68 +501,56 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, } void -WarpXParticleContainer::DepositCharge (Vector >& rho, bool local) +WarpXParticleContainer::DepositCharge (Vector >& rho, + bool local, bool reset) { int num_levels = rho.size(); int finest_level = num_levels - 1; - // each level deposits it's own particles + // each level deposits its own particles const int ng = rho[0]->nGrow(); for (int lev = 0; lev < num_levels; ++lev) { - rho[lev]->setVal(0.0, ng); - - const auto& gm = m_gdb->Geom(lev); - const auto& ba = m_gdb->ParticleBoxArray(lev); - - const Real* dx = gm.CellSize(); - const Real* plo = gm.ProbLo(); - BoxArray nba = ba; - nba.surroundingNodes(); - - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = nba[pti]; + if (reset) rho[lev]->setVal(0.0, ng); +#ifdef _OPENMP +#pragma omp parallel + { +#endif +#ifdef _OPENMP + int thread_num = omp_get_thread_num(); +#else + int thread_num = 0; +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - const auto& particles = pti.GetArrayOfStructs(); - int nstride = particles.dataShape().first; - const long np = pti.numParticles(); - FArrayBox& rhofab = (*rho[lev])[pti]; + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); + + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } - WRPX_DEPOSIT_CIC(particles.dataPtr(), nstride, np, - wp.dataPtr(), &this->charge, - rhofab.dataPtr(), box.loVect(), box.hiVect(), - plo, dx, &ng); + DepositCharge(pti, wp, ion_lev, rho[lev].get(), 0, 0, np, thread_num, lev, lev); } +#ifdef _OPENMP + } +#endif - if (!local) rho[lev]->SumBoundary(gm.periodicity()); - } +#ifdef WARPX_DIM_RZ + WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev); +#endif - // now we average down fine to crse - std::unique_ptr crse; - for (int lev = finest_level - 1; lev >= 0; --lev) { - const BoxArray& fine_BA = rho[lev+1]->boxArray(); - const DistributionMapping& fine_dm = rho[lev+1]->DistributionMap(); - BoxArray coarsened_fine_BA = fine_BA; - coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); - coarsened_fine_data.setVal(0.0); - - IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME - - for (MFIter mfi(coarsened_fine_data); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - const Box& crse_box = coarsened_fine_data[mfi].box(); - const Box& fine_box = (*rho[lev+1])[mfi].box(); - WRPX_SUM_FINE_TO_CRSE_NODAL(bx.loVect(), bx.hiVect(), ratio.getVect(), - coarsened_fine_data[mfi].dataPtr(), crse_box.loVect(), crse_box.hiVect(), - (*rho[lev+1])[mfi].dataPtr(), fine_box.loVect(), fine_box.hiVect()); + if (!local) { + const auto& gm = m_gdb->Geom(lev); + rho[lev]->SumBoundary(gm.periodicity()); } - - rho[lev]->copy(coarsened_fine_data, m_gdb->Geom(lev).periodicity(), FabArrayBase::ADD); } } @@ -841,5 +829,3 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, if (pld.m_lev == lev-1){ } } - - diff --git a/Source/WarpX.H b/Source/WarpX.H index f55670cfb..4275df3c5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -360,22 +360,6 @@ private: /// void EvolveES(int numsteps); - /// - /// Compute the gravitational potential from rho by solving Poisson's equation. - /// Both rho and phi are assumed to be node-centered. This method is only used - /// in electrostatic mode. - /// - void computePhi(const amrex::Vector >& rho, - amrex::Vector >& phi) const; - - /// - /// Compute the electric field in each direction by computing the gradient - /// the potential phi using 2nd order centered differences. Both rho and phi - /// are assumed to be node-centered. This method is only used in electrostatic mode. - /// - void computeE(amrex::Vector, 3> >& E, - const amrex::Vector >& phi) const; - // // This stuff is needed by the nodal multigrid solver when running in // electrostatic mode. @@ -411,7 +395,11 @@ private: void InitFromCheckpoint (); void PostRestart (); - void InitOpenbc (); + void InitSpaceChargeField (WarpXParticleContainer& pc); + void computePhi(const amrex::Vector >& rho, + amrex::Vector >& phi) const; + void computeE(amrex::Vector, 3> >& E, + const amrex::Vector >& phi) const; void InitPML (); void ComputePMLFactors (); -- cgit v1.2.3 From cfe63efaaf3f8e6e624d1a680655fed42d0f8236 Mon Sep 17 00:00:00 2001 From: "L. Diana Amorim" Date: Wed, 18 Dec 2019 09:55:00 -0800 Subject: Fix to RZ scalar moving window Tried to make exception message clearer --- Python/pywarpx/picmi.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) (limited to 'Python/pywarpx/picmi.py') diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 6e9c9153b..2dc08bf37 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -326,13 +326,14 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): pywarpx.geometry.prob_hi = self.upper_bound pywarpx.warpx.n_rz_azimuthal_modes = self.n_azimuthal_modes - if self.moving_window_velocity is not None and np.any(np.not_equal(self.moving_window_velocity, 0.)): - pywarpx.warpx.do_moving_window = 1 - if self.moving_window_velocity[0] != 0.: - raise Exception('In cylindrical coordinates, a moving window in r can not be done') - if self.moving_window_velocity[1] != 0.: - pywarpx.warpx.moving_window_dir = 'z' - pywarpx.warpx.moving_window_v = self.moving_window_velocity[1]/constants.c # in units of the speed of light + if self.moving_window_velocity is not None: + if np.isscalar(self.moving_window_velocity): + if self.moving_window_velocity !=0: + pywarpx.warpx.do_moving_window = 1 + pywarpx.warpx.moving_window_dir = 'z' + pywarpx.warpx.moving_window_v = self.moving_window_velocity/constants.c # in units of the speed of light + else: + raise Exception('RZ PICMI moving_window_velocity (only available in z direction) should be a scalar') if self.refined_regions: assert len(self.refined_regions) == 1, Exception('WarpX only supports one refined region.') -- cgit v1.2.3 From 32d00398e92821f093a6dd573e789d64457540df Mon Sep 17 00:00:00 2001 From: "L. Diana Amorim" Date: Wed, 18 Dec 2019 13:41:28 -0800 Subject: Corrected WarpX RZ to use PICMI moving_window_zvelocity --- Python/pywarpx/picmi.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) (limited to 'Python/pywarpx/picmi.py') diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 2dc08bf37..ce0f1f0e6 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -326,15 +326,18 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): pywarpx.geometry.prob_hi = self.upper_bound pywarpx.warpx.n_rz_azimuthal_modes = self.n_azimuthal_modes - if self.moving_window_velocity is not None: - if np.isscalar(self.moving_window_velocity): - if self.moving_window_velocity !=0: + if self.moving_window_zvelocity is not None: + if np.isscalar(self.moving_window_zvelocity): + if self.moving_window_zvelocity !=0: pywarpx.warpx.do_moving_window = 1 pywarpx.warpx.moving_window_dir = 'z' - pywarpx.warpx.moving_window_v = self.moving_window_velocity/constants.c # in units of the speed of light + pywarpx.warpx.moving_window_v = self.moving_window_zvelocity/constants.c # in units of the speed of light else: raise Exception('RZ PICMI moving_window_velocity (only available in z direction) should be a scalar') + if self.moving_window_velocity is not None: + raise Exception('PICMI RZ geometry uses moving_window_zvelocity (scalar) instead of moving_window_velocity (vector)') + if self.refined_regions: assert len(self.refined_regions) == 1, Exception('WarpX only supports one refined region.') assert self.refined_regions[0][0] == 1, Exception('The one refined region can only be level 1') -- cgit v1.2.3 From c26e7d89449586a17c9d2db1e20b916f3a484c49 Mon Sep 17 00:00:00 2001 From: "L. Diana Amorim" Date: Wed, 18 Dec 2019 13:53:10 -0800 Subject: moving_window_velocity does not exist in PICMI RZ --- Python/pywarpx/picmi.py | 3 --- 1 file changed, 3 deletions(-) (limited to 'Python/pywarpx/picmi.py') diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index ce0f1f0e6..a21775732 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -334,9 +334,6 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): pywarpx.warpx.moving_window_v = self.moving_window_zvelocity/constants.c # in units of the speed of light else: raise Exception('RZ PICMI moving_window_velocity (only available in z direction) should be a scalar') - - if self.moving_window_velocity is not None: - raise Exception('PICMI RZ geometry uses moving_window_zvelocity (scalar) instead of moving_window_velocity (vector)') if self.refined_regions: assert len(self.refined_regions) == 1, Exception('WarpX only supports one refined region.') -- cgit v1.2.3 From 7c2887709b20103d63e4a7ff0c3cfc848f17e4e5 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 20 Dec 2019 07:16:29 -0800 Subject: Fix EOL whitespace --- Python/pywarpx/picmi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Python/pywarpx/picmi.py') diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index a21775732..19d68ca29 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -334,7 +334,7 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): pywarpx.warpx.moving_window_v = self.moving_window_zvelocity/constants.c # in units of the speed of light else: raise Exception('RZ PICMI moving_window_velocity (only available in z direction) should be a scalar') - + if self.refined_regions: assert len(self.refined_regions) == 1, Exception('WarpX only supports one refined region.') assert self.refined_regions[0][0] == 1, Exception('The one refined region can only be level 1') -- cgit v1.2.3