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.

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

add depolarization field due to metal electrode to total energy density #45

Draft
wants to merge 1 commit into
base: development
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
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
13 changes: 13 additions & 0 deletions Source/FerroX.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,9 @@ AMREX_GPU_MANAGED amrex::Real FerroX::epsilonX_fe_tphase;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilonZ_fe;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_de;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_si;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_metal;
AMREX_GPU_MANAGED amrex::Real FerroX::metal_screening_length;
AMREX_GPU_MANAGED amrex::Real FerroX::metal_thickness;
AMREX_GPU_MANAGED amrex::Real FerroX::alpha; // alpha = 2*alpha_1
AMREX_GPU_MANAGED amrex::Real FerroX::beta; // beta = 4*alpha_11
AMREX_GPU_MANAGED amrex::Real FerroX::gamma; // gamma = 6*alpha_111
Expand Down Expand Up @@ -284,6 +287,16 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
pp.get("epsilonZ_fe",epsilonZ_fe);// epsilon_r for FE
pp.get("epsilon_de",epsilon_de);// epsilon_r for DE
pp.get("epsilon_si",epsilon_si);// epsilon_r for SC

epsilon_metal = 1.;
pp.query("epsilon_metal",epsilon_metal);// epsilon_r for metal

metal_screening_length = 1.e5;
pp.query("metal_screening_length",metal_screening_length);// metal_screening_length

metal_thickness = 0.;
pp.query("metal_thickness",metal_thickness);// metal_thickness

pp.get("alpha",alpha);
pp.get("beta",beta);
pp.get("gamma",FerroX::gamma);
Expand Down
3 changes: 3 additions & 0 deletions Source/FerroX_namespace.H
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@ namespace FerroX {
extern AMREX_GPU_MANAGED amrex::Real epsilonZ_fe;
extern AMREX_GPU_MANAGED amrex::Real epsilon_de;
extern AMREX_GPU_MANAGED amrex::Real epsilon_si;
extern AMREX_GPU_MANAGED amrex::Real epsilon_metal;
extern AMREX_GPU_MANAGED amrex::Real metal_screening_length;
extern AMREX_GPU_MANAGED amrex::Real metal_thickness;
extern AMREX_GPU_MANAGED amrex::Real alpha; // alpha = 2*alpha_1
extern AMREX_GPU_MANAGED amrex::Real beta; // beta = 4*alpha_11
extern AMREX_GPU_MANAGED amrex::Real gamma; // gamma = 6*alpha_111
Expand Down
7 changes: 7 additions & 0 deletions Source/Solver/TotalEnergyDensity.H
Original file line number Diff line number Diff line change
Expand Up @@ -16,3 +16,10 @@ void CalculateTDGL_RHS(Array<MultiFab, AMREX_SPACEDIM> &GL_rhs,
const Geometry& geom,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi);

void Compute_P_Sum(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum);

void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count);

void Compute_P_av(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real &P_av_z);

87 changes: 87 additions & 0 deletions Source/Solver/TotalEnergyDensity.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,23 @@ void CalculateTDGL_RHS(Array<MultiFab, AMREX_SPACEDIM> &GL_rhs,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi)
{

Real average_P_r = 0.;
Real total_P_r = 0.;
Real FE_index_counter = 0.;

Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r);

//Calculate integrated electrode charge (Qe) based on eq 13 of https://pubs.aip.org/aip/jap/article/44/8/3379/6486/Depolarization-fields-in-thin-ferroelectric-films
Real FE_thickness = FE_hi[2] - FE_lo[2];
Real coth = (exp(2.*metal_thickness/metal_screening_length) + 1.0) / (exp(2.*metal_thickness/metal_screening_length) - 1.0);
Real csch = (2.*exp(metal_thickness/metal_screening_length)) / (exp(2.*metal_thickness/metal_screening_length) - 1.0);
Real numerator = 0.5 * FE_thickness * average_P_r / epsilonX_fe;
Real denominator = metal_screening_length/epsilon_de*(coth - csch) + FE_thickness/(2.*epsilonX_fe);
Real Qe = -numerator/denominator;
Real theta = -Qe/average_P_r;
//Real E_dep_metal = -average_P_r/epsilonX_fe*(1.0 - theta);

// loop over boxes
for ( MFIter mfi(P_old[0]); mfi.isValid(); ++mfi )
{
Expand Down Expand Up @@ -153,6 +170,7 @@ void CalculateTDGL_RHS(Array<MultiFab, AMREX_SPACEDIM> &GL_rhs,
( dFdPr_Landau
+ dFdPr_grad
- Er(i,j,k)
+ pOld_r(i,j,k)/epsilonX_fe*(1.0 - theta)
//+ DphiDz(phi, z_hi, z_lo, i, j, k, dx, prob_lo, prob_hi)
);

Expand All @@ -165,4 +183,73 @@ void CalculateTDGL_RHS(Array<MultiFab, AMREX_SPACEDIM> &GL_rhs,
}
}

void Compute_P_Sum(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum)
{

// Initialize to zero
sum = 0.;

ReduceOps<ReduceOpSum> reduce_op;

ReduceData<Real> reduce_data(reduce_op);
using ReduceTuple = typename decltype(reduce_data)::Type;

for (MFIter mfi(P[2],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Box& bx_grid = mfi.validbox();

auto const& fab = P[2].array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
{
return {fab(i,j,k)};
});
}

sum = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(sum);
}


void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count)
{

// Initialize to zero
count = 0.;

ReduceOps<ReduceOpSum> reduce_op;

ReduceData<Real> reduce_data(reduce_op);
using ReduceTuple = typename decltype(reduce_data)::Type;

for (MFIter mfi(MaterialMask, TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Box& bx_grid = mfi.validbox();

auto const& fab = MaterialMask.array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
{
if(fab(i,j,k) == 0.) {
return {1.};
} else {
return {0.};
}

});
}

count = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(count);
}

void Compute_P_av(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real& P_av_z)
{
Compute_P_Sum(P, sum);
Compute_P_index_Sum(MaterialMask, count);
P_av_z = sum/count;
}