Skip to content

Commit

Permalink
add Holstrom option for handling low density regions
Browse files Browse the repository at this point in the history
  • Loading branch information
@roelof-groenewald
roelof-groenewald committed Dec 12, 2023
1 parent d97a0d3 commit 0f994a0
Show file tree
Hide file tree
Showing 4 changed files with 111 additions and 36 deletions.
9 changes: 8 additions & 1 deletion Python/pywarpx/picmi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1138,6 +1138,10 @@ class HybridPICSolver(picmistandard.base._ClassWithInit):
Whether to use the change in total current (from the first simulation
step) or the total current in the resistive term in Ohm's law.
vacuum_algorithm: str
Algorithm to use in regions where the plasma density is lower than
n_floor.
substeps: int, default=100
Number of substeps to take when updating the B-field.
Expand All @@ -1146,7 +1150,8 @@ class HybridPICSolver(picmistandard.base._ClassWithInit):
"""
def __init__(self, grid, Te=None, n0=None, gamma=None,
n_floor=None, plasma_resistivity=None,
use_dJ_in_resistive_term=None, substeps=None,
use_dJ_in_resistive_term=None,
vacuum_algorithm=None, substeps=None,
Jx_external_function=None, Jy_external_function=None,
Jz_external_function=None, **kw):
self.grid = grid
Expand All @@ -1158,6 +1163,7 @@ def __init__(self, grid, Te=None, n0=None, gamma=None,
self.n_floor = n_floor
self.plasma_resistivity = plasma_resistivity
self.use_dJ_in_resistive_term = use_dJ_in_resistive_term
self.vacuum_algorithm = vacuum_algorithm

self.substeps = substeps

Expand All @@ -1181,6 +1187,7 @@ def initialize_inputs(self):
'plasma_resistivity(rho)', self.plasma_resistivity
)
pywarpx.hybridpicmodel.use_dJ_in_resistive_term = self.use_dJ_in_resistive_term
pywarpx.hybridpicmodel.vacuum_algorithm = self.vacuum_algorithm
pywarpx.hybridpicmodel.substeps = self.substeps
pywarpx.hybridpicmodel.__setattr__(
'Jx_external_grid_function(x,y,z,t)', self.Jx_external_function
Expand Down
21 changes: 21 additions & 0 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,18 @@
#include <AMReX_REAL.H>


/**
* \brief struct to determine the algorithm used to handle low density / vacuum
regions, default is DensityFloor.
*/
struct VacuumAlgo {
enum {
DensityFloor = 0,
Holstrom = 1
};
};


/**
* \brief This class contains the parameters needed to evaluate hybrid field
* solutions (kinetic ions with fluid electrons).
Expand Down Expand Up @@ -180,6 +192,15 @@ public:
amrex::ParserExecutor<1> m_eta;
bool m_use_dJ_for_resistive_term = false;

/** Different approaches to handling low plasma density regions */
std::map<std::string, int> vacuum_algorithm_to_int = {
{"density_floor", VacuumAlgo::DensityFloor},
{"holstrom", VacuumAlgo::Holstrom},
{"default", VacuumAlgo::DensityFloor},
};
std::string m_vacuum_algorithm_str = "default";
int m_vacuum_algorithm;

/** External current */
std::string m_Jx_ext_grid_function = "0.0";
std::string m_Jy_ext_grid_function = "0.0";
Expand Down
8 changes: 8 additions & 0 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,14 @@ void HybridPICModel::ReadParameters ()
// convert electron temperature from eV to J
m_elec_temp *= PhysConst::q_e;

pp_hybrid.query("vacuum_algorithm", m_vacuum_algorithm_str);
// Convert to lower case
std::transform(
m_vacuum_algorithm_str.begin(), m_vacuum_algorithm_str.end(),
m_vacuum_algorithm_str.begin(), ::tolower
);
m_vacuum_algorithm = vacuum_algorithm_to_int[m_vacuum_algorithm_str];

// external currents
pp_hybrid.query("Jx_external_grid_function(x,y,z,t)", m_Jx_ext_grid_function);
pp_hybrid.query("Jy_external_grid_function(x,y,z,t)", m_Jy_ext_grid_function);
Expand Down
109 changes: 74 additions & 35 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -436,6 +436,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
const auto eta = hybrid_model->m_eta;
const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e;
const auto use_dJ_for_resistive_term = hybrid_model->m_use_dJ_for_resistive_term;
const auto vacuum_algorithm = hybrid_model->m_vacuum_algorithm;

// Index type required for interpolating fields from their respective
// staggering to the Ex, Ey, Ez locations
Expand Down Expand Up @@ -596,18 +597,24 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
// Skip if this cell is fully covered by embedded boundaries
if (lr(i, j, 0) <= 0) return;
#endif
// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Er_stag, coarsen, i, j, 0, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) rho_val = rho_floor;

// Get the gradient of the electron pressure
auto grad_Pe = T_Algo::UpwardDr(Pe, coefs_r, n_coefs_r, i, j, 0, 0);

// interpolate the nodal neE values to the Yee grid
auto enE_r = Interp(enE, nodal, Er_stag, coarsen, i, j, 0, 0);

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Er_stag, coarsen, i, j, 0, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) {
rho_val = rho_floor;
if (vacuum_algorithm == VacuumAlgo::Holstrom) {
enE_r = 0.0_rt;
grad_Pe = 0.0_rt;
}
}

Er(i, j, 0) = (enE_r - grad_Pe) / rho_val;

// Add resistivity only if E field value is used to update B
Expand All @@ -630,19 +637,25 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
return;
}

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Er_stag, coarsen, i, j, 0, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) rho_val = rho_floor;

// Get the gradient of the electron pressure
// -> d/dt = 0 for m = 0
auto grad_Pe = 0.0_rt;

// interpolate the nodal neE values to the Yee grid
auto enE_t = Interp(enE, nodal, Et_stag, coarsen, i, j, 0, 1);

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Et_stag, coarsen, i, j, 0, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) {
rho_val = rho_floor;
if (vacuum_algorithm == VacuumAlgo::Holstrom) {
enE_t = 0.0_rt;
grad_Pe = 0.0_rt;
}
}

Et(i, j, 0) = (enE_t - grad_Pe) / rho_val;

// Add resistivity only if E field value is used to update B
Expand All @@ -656,18 +669,24 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
// Skip field solve if this cell is fully covered by embedded boundaries
if (lz(i,j,0) <= 0) return;
#endif
// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) rho_val = rho_floor;

// Get the gradient of the electron pressure
auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, k, 0);

// interpolate the nodal neE values to the Yee grid
auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, k, 2);

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) {
rho_val = rho_floor;
if (vacuum_algorithm == VacuumAlgo::Holstrom) {
enE_z = 0.0_rt;
grad_Pe = 0.0_rt;
}
}

Ez(i, j, k) = (enE_z - grad_Pe) / rho_val;

// Add resistivity only if E field value is used to update B
Expand Down Expand Up @@ -714,6 +733,7 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
const auto eta = hybrid_model->m_eta;
const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e;
const auto use_dJ_for_resistive_term = hybrid_model->m_use_dJ_for_resistive_term;
const auto vacuum_algorithm = hybrid_model->m_vacuum_algorithm;

// Index type required for interpolating fields from their respective
// staggering to the Ex, Ey, Ez locations
Expand Down Expand Up @@ -872,18 +892,24 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
// Skip if this cell is fully covered by embedded boundaries
if (lx(i, j, k) <= 0) return;
#endif
// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ex_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) rho_val = rho_floor;

// Get the gradient of the electron pressure
auto grad_Pe = T_Algo::UpwardDx(Pe, coefs_x, n_coefs_x, i, j, k);

// interpolate the nodal neE values to the Yee grid
auto enE_x = Interp(enE, nodal, Ex_stag, coarsen, i, j, k, 0);

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ex_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) {
rho_val = rho_floor;
if (vacuum_algorithm == VacuumAlgo::Holstrom) {
enE_x = 0.0_rt;
grad_Pe = 0.0_rt;
}
}

Ex(i, j, k) = (enE_x - grad_Pe) / rho_val;

// Add resistivity only if E field value is used to update B
Expand All @@ -903,18 +929,24 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
if (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j-1, k)<=0 || lz(i, j, k)<=0) return;
#endif
#endif
// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ey_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) rho_val = rho_floor;

// Get the gradient of the electron pressure
auto grad_Pe = T_Algo::UpwardDy(Pe, coefs_y, n_coefs_y, i, j, k);

// interpolate the nodal neE values to the Yee grid
auto enE_y = Interp(enE, nodal, Ey_stag, coarsen, i, j, k, 1);

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ey_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) {
rho_val = rho_floor;
if (vacuum_algorithm == VacuumAlgo::Holstrom) {
enE_y = 0.0_rt;
grad_Pe = 0.0_rt;
}
}

Ey(i, j, k) = (enE_y - grad_Pe) / rho_val;

// Add resistivity only if E field value is used to update B
Expand All @@ -928,18 +960,25 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
// Skip field solve if this cell is fully covered by embedded boundaries
if (lz(i,j,k) <= 0) return;
#endif
// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) rho_val = rho_floor;

// Get the gradient of the electron pressure
auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, k);

// interpolate the nodal neE values to the Yee grid
auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, k, 2);

// Interpolate to get the appropriate charge density in space
Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0);

// safety condition since we divide by rho_val later
if (rho_val < rho_floor) {
rho_val = rho_floor;
if (vacuum_algorithm == VacuumAlgo::Holstrom) {
enE_z = 0.0_rt;
grad_Pe = 0.0_rt;
}
}

Ez(i, j, k) = (enE_z - grad_Pe) / rho_val;

// Add resistivity only if E field value is used to update B
Expand Down

0 comments on commit 0f994a0

Please sign in to comment.