Skip to content

Added JJ function #152

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

Draft
tyh0123 wants to merge 1 commit into
base: development
Choose a base branch
from
+579 −30
Draft

Added JJ function #152

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion GNUmakefile
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
AMREX_HOME ?= ../amrex
AMREX_HOME = ../../amrex

DEBUG = FALSE
#DEBUG = TRUE
Expand Down
12 changes: 12 additions & 0 deletions Source/Diagnostics/FullDiagnostics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -716,6 +716,18 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getInductor().m_inductor_y_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "inductorz" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getInductor().m_inductor_z_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "josephson_phi_x" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getJosephsonJunction().m_phi_x_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "josephson_phi_y" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getJosephsonJunction().m_phi_y_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "josephson_phi_z" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getJosephsonJunction().m_phi_z_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "josephson_Ic_x" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getJosephsonJunction().m_Ic_x_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "josephson_Ic_y" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getJosephsonJunction().m_Ic_y_mf.get(), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "josephson_Ic_z" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getJosephsonJunction().m_Ic_z_mf.get(), lev, m_crse_ratio);
}
else {

Expand Down
8 changes: 8 additions & 0 deletions Source/Evolve/WarpXEvolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,6 +154,10 @@ WarpX::Evolve (int numsteps)
m_inductor->EvolveInductorJ(-0.5_rt*dt[0]); // J^(n) to J^(n-1/2) using E^(n)
FillBoundaryJ(guard_cells.ng_alloc_EB);
}
if (use_josephson_junction) {
m_jj->EvolveJunctionJ(-0.5_rt*dt[0]); // phi^(0) to phi^(-1/2), add JJ J^(-1/2)
FillBoundaryJ(guard_cells.ng_alloc_EB);
}
is_synchronized = false;
} else {
if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) {
Expand Down Expand Up @@ -431,6 +435,10 @@ WarpX::OneStep_nosub (Real cur_time)
m_inductor->EvolveInductorJ(dt[0]); // J^(n-1/2) to J^(n+1/2) using E^(n)
FillBoundaryJ(guard_cells.ng_alloc_EB);
}
if (use_josephson_junction == 1) {
m_jj->EvolveJunctionJ(dt[0]); // phi^(n-1/2) to phi^(n+1/2), add JJ J^(n+1/2)
FillBoundaryJ(guard_cells.ng_alloc_EB);
}

ExecutePythonCallback("afterdeposition");

Expand Down
104 changes: 80 additions & 24 deletions Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,12 +133,37 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
amrex::ignore_unused(edge_lengths);
#endif

auto & warpx = WarpX::GetInstance();
const int lev = 0;
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = warpx.Geom(lev).CellSizeArray();

const int use_lumped_resistor = warpx.use_lumped_resistor;
const int use_lumped_capacitor = warpx.use_lumped_capacitor;

amrex::MultiFab& sigma_mf = macroscopic_properties->getsigma_mf();
amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf();
#ifndef WARPX_MAG_LLG
amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf();
#endif

amrex::MultiFab* lumped_resistor_x_mf = nullptr;
amrex::MultiFab* lumped_resistor_y_mf = nullptr;
amrex::MultiFab* lumped_resistor_z_mf = nullptr;
amrex::MultiFab* lumped_capacitor_x_mf = nullptr;
amrex::MultiFab* lumped_capacitor_y_mf = nullptr;
amrex::MultiFab* lumped_capacitor_z_mf = nullptr;

if (use_lumped_resistor) {
lumped_resistor_x_mf = macroscopic_properties->get_pointer_lumped_resistor_x();
lumped_resistor_y_mf = macroscopic_properties->get_pointer_lumped_resistor_y();
lumped_resistor_z_mf = macroscopic_properties->get_pointer_lumped_resistor_z();
}
if (use_lumped_capacitor) {
lumped_capacitor_x_mf = macroscopic_properties->get_pointer_lumped_capacitor_x();
lumped_capacitor_y_mf = macroscopic_properties->get_pointer_lumped_capacitor_y();
lumped_capacitor_z_mf = macroscopic_properties->get_pointer_lumped_capacitor_z();
}

// Index type required for calling ablastr::coarsen::sample::Interp to interpolate macroscopic
// properties from their respective staggering to the Ex, Ey, Ez locations
amrex::GpuArray<int, 3> const& sigma_stag = macroscopic_properties->sigma_IndexType;
Expand Down Expand Up @@ -180,6 +205,20 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
amrex::Array4<amrex::Real> const& mu_arr = mu_mf.array(mfi);
#endif

// Lumped element arrays (default-constructed; populated only when feature enabled).
amrex::Array4<amrex::Real> resistor_x_arr, resistor_y_arr, resistor_z_arr;
amrex::Array4<amrex::Real> capacitor_x_arr, capacitor_y_arr, capacitor_z_arr;
if (use_lumped_resistor) {
resistor_x_arr = lumped_resistor_x_mf->array(mfi);
resistor_y_arr = lumped_resistor_y_mf->array(mfi);
resistor_z_arr = lumped_resistor_z_mf->array(mfi);
}
if (use_lumped_capacitor) {
capacitor_x_arr = lumped_capacitor_x_mf->array(mfi);
capacitor_y_arr = lumped_capacitor_y_mf->array(mfi);
capacitor_z_arr = lumped_capacitor_z_mf->array(mfi);
}

// Extract stencil coefficients
Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
int const n_coefs_x = m_stencil_coefs_x.size();
Expand Down Expand Up @@ -214,18 +253,25 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
// Skip field push if this cell is fully covered by embedded boundaries
if (lx(i, j, k) <= 0) return;
#endif
// Interpolate conductivity, sigma, to Ex position on the grid
amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag,
Ex_stag, macro_cr, i, j, k, scomp);
// Interpolated permittivity, epsilon, to Ex position on the grid
amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag,
Ex_stag, macro_cr, i, j, k, scomp);
amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt);
amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt);
Ex(i, j, k) = alpha * Ex(i, j, k)
+ beta * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0)
+ T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0)
) - beta * jx(i, j, k);
// fac1 = (dt/eps)*(sigma + lumped-resistor equivalent conductance)
amrex::Real const fac1 = (dt/epsilon_interp) *
( sigma_interp + ((use_lumped_resistor && resistor_x_arr(i,j,k)!=0._rt)
? dx[0] / (dx[1]*dx[2]*resistor_x_arr(i,j,k))
: 0._rt) );
// fac2 = lumped-capacitor / natural cell capacitance
amrex::Real const fac2 = (use_lumped_capacitor && capacitor_x_arr(i,j,k)!=0._rt)
? capacitor_x_arr(i,j,k) * dx[0] / (dx[1]*dx[2]*epsilon_interp)
: 0._rt;
amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2);
amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2);
Ex(i, j, k) = alpha_c * Ex(i, j, k)
+ (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0)
+ T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0))
- (dt/epsilon_interp) * beta_c * jx(i, j, k);
},

[=] AMREX_GPU_DEVICE (int i, int j, int k){
Expand All @@ -238,37 +284,47 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
if (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j, k)<=0 || lz(i, j-1, k)<=0) return;
#endif
#endif
// Interpolate conductivity, sigma, to Ey position on the grid
amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag,
Ey_stag, macro_cr, i, j, k, scomp);
// Interpolated permittivity, epsilon, to Ey position on the grid
amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag,
Ey_stag, macro_cr, i, j, k, scomp);
amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt);
amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt);
Ey(i, j, k) = alpha * Ey(i, j, k)
+ beta * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0)
+ T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0)
) - beta * jy(i, j, k);
amrex::Real const fac1 = (dt/epsilon_interp) *
( sigma_interp + ((use_lumped_resistor && resistor_y_arr(i,j,k)!=0._rt)
? dx[1] / (dx[0]*dx[2]*resistor_y_arr(i,j,k))
: 0._rt) );
amrex::Real const fac2 = (use_lumped_capacitor && capacitor_y_arr(i,j,k)!=0._rt)
? capacitor_y_arr(i,j,k) * dx[1] / (dx[0]*dx[2]*epsilon_interp)
: 0._rt;
amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2);
amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2);
Ey(i, j, k) = alpha_c * Ey(i, j, k)
+ (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0)
+ T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0))
- (dt/epsilon_interp) * beta_c * jy(i, j, k);
},

[=] AMREX_GPU_DEVICE (int i, int j, int k){
#ifdef AMREX_USE_EB
// Skip field push if this cell is fully covered by embedded boundaries
if (lz(i,j,k) <= 0) return;
#endif
// Interpolate conductivity, sigma, to Ez position on the grid
amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag,
Ez_stag, macro_cr, i, j, k, scomp);
// Interpolated permittivity, epsilon, to Ez position on the grid
amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag,
Ez_stag, macro_cr, i, j, k, scomp);
amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt);
amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt);
Ez(i, j, k) = alpha * Ez(i, j, k)
+ beta * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0)
+ T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0)
) - beta * jz(i, j, k);
amrex::Real const fac1 = (dt/epsilon_interp) *
( sigma_interp + ((use_lumped_resistor && resistor_z_arr(i,j,k)!=0._rt)
? dx[2] / (dx[0]*dx[1]*resistor_z_arr(i,j,k))
: 0._rt) );
amrex::Real const fac2 = (use_lumped_capacitor && capacitor_z_arr(i,j,k)!=0._rt)
? capacitor_z_arr(i,j,k) * dx[2] / (dx[0]*dx[1]*epsilon_interp)
: 0._rt;
amrex::Real const alpha_c = T_MacroAlgo::alpha_compact(fac1, fac2);
amrex::Real const beta_c = T_MacroAlgo::beta_compact(fac1, fac2);
Ez(i, j, k) = alpha_c * Ez(i, j, k)
+ (dt/epsilon_interp) * beta_c * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0)
+ T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0))
- (dt/epsilon_interp) * beta_c * jz(i, j, k);
}
);
}
Expand Down
91 changes: 91 additions & 0 deletions Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,22 @@ public:
amrex::MultiFab& getmu_mf () {return (*m_mu_mf);}
amrex::MultiFab * get_pointer_mu () {return m_mu_mf.get();}

/** return MultiFabs holding the per-edge lumped resistor values [Ohm] */
amrex::MultiFab& getlumped_resistor_x_mf () {return (*m_lumped_resistor_x_mf);}
amrex::MultiFab * get_pointer_lumped_resistor_x () {return m_lumped_resistor_x_mf.get();}
amrex::MultiFab& getlumped_resistor_y_mf () {return (*m_lumped_resistor_y_mf);}
amrex::MultiFab * get_pointer_lumped_resistor_y () {return m_lumped_resistor_y_mf.get();}
amrex::MultiFab& getlumped_resistor_z_mf () {return (*m_lumped_resistor_z_mf);}
amrex::MultiFab * get_pointer_lumped_resistor_z () {return m_lumped_resistor_z_mf.get();}

/** return MultiFabs holding the per-edge lumped capacitor values [Farad] */
amrex::MultiFab& getlumped_capacitor_x_mf () {return (*m_lumped_capacitor_x_mf);}
amrex::MultiFab * get_pointer_lumped_capacitor_x () {return m_lumped_capacitor_x_mf.get();}
amrex::MultiFab& getlumped_capacitor_y_mf () {return (*m_lumped_capacitor_y_mf);}
amrex::MultiFab * get_pointer_lumped_capacitor_y () {return m_lumped_capacitor_y_mf.get();}
amrex::MultiFab& getlumped_capacitor_z_mf () {return (*m_lumped_capacitor_z_mf);}
amrex::MultiFab * get_pointer_lumped_capacitor_z () {return m_lumped_capacitor_z_mf.get();}

/** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */
amrex::GpuArray<int, 3> macro_cr_ratio;
/** Initializes the Multifabs storing macroscopic properties
Expand Down Expand Up @@ -86,6 +102,21 @@ public:
/** Gpu Vector with index type of the Bz multifab */
amrex::GpuArray<int, 3> Bz_IndexType;

/** Gpu Vector with index type of the lumped resistor multifab */
amrex::GpuArray<int, 3> lumped_resistor_x_IndexType;
amrex::GpuArray<int, 3> lumped_resistor_y_IndexType;
amrex::GpuArray<int, 3> lumped_resistor_z_IndexType;

/** Gpu Vector with index type of the lumped capacitor multifab */
amrex::GpuArray<int, 3> lumped_capacitor_x_IndexType;
amrex::GpuArray<int, 3> lumped_capacitor_y_IndexType;
amrex::GpuArray<int, 3> lumped_capacitor_z_IndexType;

/** Gpu Vector with index type of jx/jy/jz multifab (edge-centered) */
amrex::GpuArray<int, 3> jx_IndexType;
amrex::GpuArray<int, 3> jy_IndexType;
amrex::GpuArray<int, 3> jz_IndexType;

/** Stores initialization type for conductivity : constant or parser */
std::string m_sigma_s = "constant";
/** Stores initialization type for permittivity : constant or parser */
Expand Down Expand Up @@ -117,6 +148,14 @@ public:
std::unique_ptr<amrex::Parser> m_epsilon_parser;
std::unique_ptr<amrex::Parser> m_mu_parser;

std::unique_ptr<amrex::Parser> m_lumped_resistor_x_parser;
std::unique_ptr<amrex::Parser> m_lumped_resistor_y_parser;
std::unique_ptr<amrex::Parser> m_lumped_resistor_z_parser;

std::unique_ptr<amrex::Parser> m_lumped_capacitor_x_parser;
std::unique_ptr<amrex::Parser> m_lumped_capacitor_y_parser;
std::unique_ptr<amrex::Parser> m_lumped_capacitor_z_parser;

#ifdef WARPX_MAG_LLG
/** Gpu Vector with index type of the Hx multifab */
amrex::GpuArray<int, 3> Hx_IndexType;
Expand Down Expand Up @@ -330,6 +369,22 @@ private:
std::string m_str_epsilon_function;
std::string m_str_mu_function;

/** Multifabs for lumped resistor (per-edge Ohm values) */
std::unique_ptr<amrex::MultiFab> m_lumped_resistor_x_mf;
std::unique_ptr<amrex::MultiFab> m_lumped_resistor_y_mf;
std::unique_ptr<amrex::MultiFab> m_lumped_resistor_z_mf;
std::string m_str_lumped_resistor_x_function;
std::string m_str_lumped_resistor_y_function;
std::string m_str_lumped_resistor_z_function;

/** Multifabs for lumped capacitor (per-edge Farad values) */
std::unique_ptr<amrex::MultiFab> m_lumped_capacitor_x_mf;
std::unique_ptr<amrex::MultiFab> m_lumped_capacitor_y_mf;
std::unique_ptr<amrex::MultiFab> m_lumped_capacitor_z_mf;
std::string m_str_lumped_capacitor_x_function;
std::string m_str_lumped_capacitor_y_function;
std::string m_str_lumped_capacitor_z_function;

};

/**
Expand Down Expand Up @@ -362,6 +417,24 @@ struct LaxWendroffAlgo {
return beta;
}

/** Compact form for combined dissipation (fac1) and added capacitance (fac2).
* Used when lumped resistor/capacitor terms are present. */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real alpha_compact (amrex::Real const fac1,
amrex::Real const fac2) {
using namespace amrex;
amrex::Real alpha_compact = (1._rt - 0.5_rt * fac1 + fac2)/(1._rt + 0.5_rt * fac1 + fac2);
return alpha_compact;
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real beta_compact (amrex::Real const fac1,
amrex::Real const fac2) {
using namespace amrex;
amrex::Real beta_compact = 1._rt / (1._rt + 0.5_rt * fac1 + fac2);
return beta_compact;
}

};

/**
Expand Down Expand Up @@ -394,6 +467,24 @@ struct BackwardEulerAlgo {
return beta;
}

/** Compact form for combined dissipation (fac1) and added capacitance (fac2).
* Used when lumped resistor/capacitor terms are present. */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real alpha_compact (amrex::Real const fac1,
amrex::Real const fac2) {
using namespace amrex;
amrex::Real alpha_compact = (1._rt + fac2)/(1._rt + fac1 + fac2);
return alpha_compact;
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real beta_compact (amrex::Real const fac1,
amrex::Real const fac2) {
using namespace amrex;
amrex::Real beta_compact = 1._rt / (1._rt + fac1 + fac2);
return beta_compact;
}

};

#endif // WARPX_MACROSCOPIC_PROPERTIES_H_
Loading