From ce30a1f93500659b3f0f05390b9b3f45ec0e749c Mon Sep 17 00:00:00 2001
From: Prabhat Kumar <prabhatkumar@lbl.gov>
Date: Thu, 2 Nov 2023 15:22:07 -0700
Subject: [PATCH] add depolarization field due to metal electrode to total
 energy density

---
 Source/FerroX.cpp                    | 13 +++++
 Source/FerroX_namespace.H            |  3 +
 Source/Solver/TotalEnergyDensity.H   |  7 +++
 Source/Solver/TotalEnergyDensity.cpp | 87 ++++++++++++++++++++++++++++
 4 files changed, 110 insertions(+)

diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp
index 680074f..53c31bf 100644
--- a/Source/FerroX.cpp
+++ b/Source/FerroX.cpp
@@ -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
@@ -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);
diff --git a/Source/FerroX_namespace.H b/Source/FerroX_namespace.H
index d5d60fb..0c299e1 100644
--- a/Source/FerroX_namespace.H
+++ b/Source/FerroX_namespace.H
@@ -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
diff --git a/Source/Solver/TotalEnergyDensity.H b/Source/Solver/TotalEnergyDensity.H
index 15a1326..e15c0a9 100644
--- a/Source/Solver/TotalEnergyDensity.H
+++ b/Source/Solver/TotalEnergyDensity.H
@@ -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);
+
diff --git a/Source/Solver/TotalEnergyDensity.cpp b/Source/Solver/TotalEnergyDensity.cpp
index ce2cb01..df57dc5 100644
--- a/Source/Solver/TotalEnergyDensity.cpp
+++ b/Source/Solver/TotalEnergyDensity.cpp
@@ -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 )
         {
@@ -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)
                     );
 
@@ -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;
+}