Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[WIP] Avoid side-effect: Undo the scaling of rho in the computePhi #3445

Open
wants to merge 2 commits into
base: development
Choose a base branch
from
Open

[WIP] Avoid side-effect: Undo the scaling of rho in the computePhi #3445

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b2e5a33
Undo the scaling of rho in the `computePhi`
RemiLehe Oct 3, 2022
32e6816
Modify sigma to include -epsilon0
RemiLehe Oct 3, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 28 additions & 6 deletions Source/ablastr/fields/PoissonSolver.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,6 @@ computePhi (amrex::Vector<amrex::MultiFab*> const & rho,
// scale rho appropriately; also determine if rho is zero everywhere
amrex::Real max_norm_b = 0.0;
for (int lev=0; lev<=finest_level; lev++) {
rho[lev]->mult(-1._rt/PhysConst::ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect!
max_norm_b = amrex::max(max_norm_b, rho[lev]->norm0());
}
amrex::ParallelDescriptor::ReduceRealMax(max_norm_b);
Expand Down Expand Up @@ -226,12 +225,14 @@ computePhi (amrex::Vector<amrex::MultiFab*> const & rho,
// one of the axes of the grid, i.e. that only *one* of the
// components of `beta` is non-negligible. // we use this
#if defined(WARPX_DIM_RZ)
linop.setSigma({0._rt, 1._rt-beta_solver[1]*beta_solver[1]});
linop.setSigma({
-PhysConst::ep0,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you're intending to multiply through by ep0, this looks like a mistake - the first entry would still be 0, wouldn't it?

-PhysConst::ep0 * (1._rt-beta_solver[1]*beta_solver[1])});
#else
linop.setSigma({AMREX_D_DECL(
1._rt-beta_solver[0]*beta_solver[0],
1._rt-beta_solver[1]*beta_solver[1],
1._rt-beta_solver[2]*beta_solver[2])});
-PhysConst::ep0 * (1._rt-beta_solver[0]*beta_solver[0]),
-PhysConst::ep0 * (1._rt-beta_solver[1]*beta_solver[1]),
-PhysConst::ep0 * (1._rt-beta_solver[2]*beta_solver[2]))});
#endif

#if defined(AMREX_USE_EB)
Expand All @@ -243,12 +244,33 @@ computePhi (amrex::Vector<amrex::MultiFab*> const & rho,
else
linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value()));
#endif

#else
// In the absence of EB and RZ: use a more generic solver
// that can handle beams propagating in any direction
amrex::MLNodeTensorLaplacian linop( {geom[lev]}, {grids[lev]},
{dmap[lev]}, info );
linop.setBeta( beta_solver ); // for the non-axis-aligned solver

#if defined(WARPX_DIM_1D_Z)
linop.setSigma({
-PhysConst::ep0 * (1._rt-beta_solver[0]*beta_solver[0]), //xx
});
#elif defined(WARPX_DIM_XZ)
linop.setSigma({
-PhysConst::ep0 * (1._rt-beta_solver[0]*beta_solver[0]), //xx
-PhysConst::ep0 * (-beta_solver[0]*beta_solver[1]), //xy
-PhysConst::ep0 * (1._rt-beta_solver[1]*beta_solver[1]) //yy
});
#else
linop.setSigma({
-PhysConst::ep0 * (1._rt-beta_solver[0]*beta_solver[0]), //xx
-PhysConst::ep0 * (-beta_solver[0]*beta_solver[1]), //xy
-PhysConst::ep0 * (-beta_solver[0]*beta_solver[2]), //xz
-PhysConst::ep0 * (1._rt-beta_solver[1]*beta_solver[1]), //yy
-PhysConst::ep0 * (1._rt-beta_solver[1]*beta_solver[2]), //yz
-PhysConst::ep0 * (1._rt-beta_solver[2]*beta_solver[2]) //zz
});
#endif
#endif

// Solve the Poisson equation
Expand Down