From feade4fd966b6eb032cab695c8e45476236b796e Mon Sep 17 00:00:00 2001 From: prkkumar Date: Wed, 2 Aug 2023 09:19:22 -0700 Subject: [PATCH 1/5] drift-diffusion model for carrier transport --- Source/FerroX.cpp | 16 +++ Source/FerroX_namespace.H | 4 + Source/Solver/ChargeDensity.H | 8 ++ Source/Solver/ChargeDensity.cpp | 171 ++++++++++++++++++++++++++++ Source/Solver/DerivativeAlgorithm.H | 9 ++ Source/main.cpp | 19 +++- 6 files changed, 224 insertions(+), 3 deletions(-) diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 861ba14..e8fcb60 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -207,6 +207,12 @@ AMREX_GPU_MANAGED amrex::Real FerroX::T; AMREX_GPU_MANAGED amrex::Real FerroX::acceptor_doping; AMREX_GPU_MANAGED amrex::Real FerroX::donor_doping; +//Drift-Diffusion +AMREX_GPU_MANAGED amrex::Real FerroX::electron_mobility; +AMREX_GPU_MANAGED amrex::Real FerroX::hole_mobility; +AMREX_GPU_MANAGED amrex::Real FerroX::electron_affinity; +AMREX_GPU_MANAGED amrex::Real FerroX::bandgap; + // P and Phi Bc AMREX_GPU_MANAGED amrex::Real FerroX::lambda; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_lo; @@ -440,6 +446,16 @@ void InitializeFerroXNamespace(const amrex::GpuArray &Jn, + Array &Jp, + MultiFab& rho, + MultiFab& e_den, + MultiFab& p_den, + const Geometry& geom, + const MultiFab& MaterialMask); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index f418a6f..fc013ad 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -1,4 +1,5 @@ #include "ChargeDensity.H" +#include "DerivativeAlgorithm.H" // Compute rho in SC region for given phi void ComputeRho(MultiFab& PoissonPhi, @@ -70,3 +71,173 @@ void ComputeRho(MultiFab& PoissonPhi, } } +// Drift-Diffusion : Compute rho in SC region for given phi, Jn, and Jp +void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, + Array &Jn, + Array &Jp, + MultiFab& rho, + MultiFab& e_den, + MultiFab& p_den, + const Geometry& geom, + const MultiFab& MaterialMask) +{ + MultiFab Ec_mf(rho.boxArray(), rho.DistributionMap(), 1, 0); + MultiFab Ev_mf(rho.boxArray(), rho.DistributionMap(), 1, 0); + + //Calculate Ec and Ev + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + const Array4& phi = PoissonPhi.array(mfi); + const Array4& Ec_arr = Ec_mf.array(mfi); + const Array4& Ev_arr = Ev_mf.array(mfi); + const Array4& mask = MaterialMask.array(mfi); + + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + + if (mask(i,j,k) >= 2.0) { + //electrons + Ec_arr(i,j,k) = -q*phi(i,j,k) - electron_affinity; + Ev_arr(i,j,k) = Ec_arr(i,j,k) - bandgap; + } else { + Ec_arr(i,j,k) = 0.0; + Ev_arr(i,j,k) = 0.0; + } + }); + } + + // Calculate currents Jn and Jp + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); + + // Calculate charge density from Phi, Nc, Nv, Ec, and Ev + + const Array4& hole_den_arr = p_den.array(mfi); + const Array4& e_den_arr = e_den.array(mfi); + const Array4& phi = PoissonPhi.array(mfi); + const Array4& mask = MaterialMask.array(mfi); + + const Array4& Ec_arr = Ec_mf.array(mfi); + const Array4& Ev_arr = Ev_mf.array(mfi); + const Array4& Jnx_arr = Jn[0].array(mfi); + const Array4& Jny_arr = Jn[1].array(mfi); + const Array4& Jnz_arr = Jn[2].array(mfi); + + const Array4& Jpx_arr = Jp[0].array(mfi); + const Array4& Jpy_arr = Jp[1].array(mfi); + const Array4& Jpz_arr = Jp[2].array(mfi); + + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + + if (mask(i,j,k) >= 2.0) { + + //electrons + Real gradEc_x = DFDx(Ec_arr, i, j, k, dx); + Real gradEc_y = DFDy(Ec_arr, i, j, k, dx); + Real gradEc_z = DFDz(Ec_arr, i, j, k, dx); + + Real gradn_x = DFDx(e_den_arr, i, j, k, dx); + Real gradn_y = DFDy(e_den_arr, i, j, k, dx); + Real gradn_z = DFDz(e_den_arr, i, j, k, dx); + + Jnx_arr(i,j,k) = electron_mobility*e_den_arr(i,j,k)*gradEc_x + kb*T*electron_mobility*gradn_x; + Jny_arr(i,j,k) = electron_mobility*e_den_arr(i,j,k)*gradEc_y + kb*T*electron_mobility*gradn_y; + Jnz_arr(i,j,k) = electron_mobility*e_den_arr(i,j,k)*gradEc_z + kb*T*electron_mobility*gradn_z; + + //holes + Real gradEv_x = DFDx(Ev_arr, i, j, k, dx); + Real gradEv_y = DFDy(Ev_arr, i, j, k, dx); + Real gradEv_z = DFDz(Ev_arr, i, j, k, dx); + + Real gradp_x = DFDx(hole_den_arr, i, j, k, dx); + Real gradp_y = DFDy(hole_den_arr, i, j, k, dx); + Real gradp_z = DFDz(hole_den_arr, i, j, k, dx); + + Jpx_arr(i,j,k) = hole_mobility*hole_den_arr(i,j,k)*gradEv_x - kb*T*hole_mobility*gradp_x; + Jpy_arr(i,j,k) = hole_mobility*hole_den_arr(i,j,k)*gradEv_y - kb*T*hole_mobility*gradp_y; + Jpz_arr(i,j,k) = hole_mobility*hole_den_arr(i,j,k)*gradEv_z - kb*T*hole_mobility*gradp_z; + + } else { + + Jnx_arr(i,j,k) = 0.0; + Jny_arr(i,j,k) = 0.0; + Jnz_arr(i,j,k) = 0.0; + + Jpx_arr(i,j,k) = 0.0; + Jpy_arr(i,j,k) = 0.0; + Jpz_arr(i,j,k) = 0.0; + + } + }); + } + + // loop over boxes + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); + + // Calculate charge density from Phi, Nc, Nv, Ec, and Ev + MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); + MultiFab donor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); + MultiFab div_Jn(rho.boxArray(), rho.DistributionMap(), 1, 0); + MultiFab div_Jp(rho.boxArray(), rho.DistributionMap(), 1, 0); + + const Array4& hole_den_arr = p_den.array(mfi); + const Array4& e_den_arr = e_den.array(mfi); + const Array4& charge_den_arr = rho.array(mfi); + const Array4& phi = PoissonPhi.array(mfi); + const Array4& acceptor_den_arr = acceptor_den.array(mfi); + const Array4& donor_den_arr = donor_den.array(mfi); + const Array4& mask = MaterialMask.array(mfi); + + const Array4& Jnx_arr = Jn[0].array(mfi); + const Array4& Jny_arr = Jn[1].array(mfi); + const Array4& Jnz_arr = Jn[2].array(mfi); + const Array4& div_Jn_arr = div_Jn.array(mfi); + + const Array4& Jpx_arr = Jp[0].array(mfi); + const Array4& Jpy_arr = Jp[1].array(mfi); + const Array4& Jpz_arr = Jp[2].array(mfi); + const Array4& div_Jp_arr = div_Jp.array(mfi); + + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + + if (mask(i,j,k) >= 2.0) { + + div_Jn_arr(i,j,k) = DFDx(Jnx_arr, i,j,k,dx) + DFDy(Jny_arr, i,j,k,dx) + DFDz(Jnz_arr, i,j,k,dx); + div_Jp_arr(i,j,k) = DFDx(Jpx_arr, i,j,k,dx) + DFDy(Jpy_arr, i,j,k,dx) + DFDz(Jpz_arr, i,j,k,dx); + + e_den_arr(i,j,k) += dt*div_Jn_arr(i,j,k); + hole_den_arr(i,j,k) -= dt*div_Jp_arr(i,j,k); + + //If in channel, set acceptor doping, else (Source/Drain) set donor doping + if (mask(i,j,k) == 3.0) { + acceptor_den_arr(i,j,k) = acceptor_doping; + donor_den_arr(i,j,k) = 0.0; + } else { // Source / Drain + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = donor_doping; + } + + charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + + } else { + + charge_den_arr(i,j,k) = 0.0; + + } + }); + } + } + diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H index 8f3256d..f01fc93 100644 --- a/Source/Solver/DerivativeAlgorithm.H +++ b/Source/Solver/DerivativeAlgorithm.H @@ -39,6 +39,15 @@ using namespace FerroX; return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]); } + /** + * Perform first derivative dF/dz */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DFDz ( + amrex::Array4 const& F, + int const i, int const j, int const k, amrex::GpuArray dx) { + return (F(i,j,k+1) - F(i,j,k-1))/(2.*dx[2]); + } + /** * Perform first derivative dP/dx */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/main.cpp b/Source/main.cpp index a9e6740..bcf567a 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -125,6 +125,15 @@ void main_main (c_FerroX& rFerroX) MultiFab angle_beta(ba, dm, 1, 0); MultiFab angle_theta(ba, dm, 1, 0); + //Drift-Diffusion + Array Jn; + Array Jp; + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) + { + Jn[dir].define(ba, dm, Ncomp, 0); + Jp[dir].define(ba, dm, Ncomp, 0); + } + //Initialize material mask InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); //InitializeMaterialMask(rFerroX, geom, MaterialMask); @@ -403,7 +412,9 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -473,7 +484,8 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -609,7 +621,8 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates From be918146742a15b0b0a64de74e6d482ffb29292c Mon Sep 17 00:00:00 2001 From: prkkumar Date: Wed, 2 Aug 2023 12:45:07 -0700 Subject: [PATCH 2/5] fix time integration of continuity equation --- Source/Solver/ChargeDensity.H | 2 ++ Source/Solver/ChargeDensity.cpp | 8 ++++-- Source/main.cpp | 44 ++++++++++++++++++++++++--------- 3 files changed, 40 insertions(+), 14 deletions(-) diff --git a/Source/Solver/ChargeDensity.H b/Source/Solver/ChargeDensity.H index 1756324..0ef978f 100644 --- a/Source/Solver/ChargeDensity.H +++ b/Source/Solver/ChargeDensity.H @@ -18,5 +18,7 @@ void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, const Geometry& geom, const MultiFab& MaterialMask); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index fc013ad..dc9a261 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -78,6 +78,8 @@ void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, const Geometry& geom, const MultiFab& MaterialMask) { @@ -194,6 +196,8 @@ void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, const Array4& hole_den_arr = p_den.array(mfi); const Array4& e_den_arr = e_den.array(mfi); + const Array4& hole_den_old_arr = p_den.array(mfi); + const Array4& e_den_old_arr = e_den.array(mfi); const Array4& charge_den_arr = rho.array(mfi); const Array4& phi = PoissonPhi.array(mfi); const Array4& acceptor_den_arr = acceptor_den.array(mfi); @@ -218,8 +222,8 @@ void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, div_Jn_arr(i,j,k) = DFDx(Jnx_arr, i,j,k,dx) + DFDy(Jny_arr, i,j,k,dx) + DFDz(Jnz_arr, i,j,k,dx); div_Jp_arr(i,j,k) = DFDx(Jpx_arr, i,j,k,dx) + DFDy(Jpy_arr, i,j,k,dx) + DFDz(Jpz_arr, i,j,k,dx); - e_den_arr(i,j,k) += dt*div_Jn_arr(i,j,k); - hole_den_arr(i,j,k) -= dt*div_Jp_arr(i,j,k); + e_den_arr(i,j,k) = e_den_old_arr(i,j,k) + dt*div_Jn_arr(i,j,k); + hole_den_arr(i,j,k) = hole_den_old_arr(i,j,k) - dt*div_Jp_arr(i,j,k); //If in channel, set acceptor doping, else (Source/Drain) set donor doping if (mask(i,j,k) == 3.0) { diff --git a/Source/main.cpp b/Source/main.cpp index bcf567a..dacd46c 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -118,6 +118,8 @@ void main_main (c_FerroX& rFerroX) MultiFab hole_den(ba, dm, 1, 0); MultiFab e_den(ba, dm, 1, 0); + MultiFab hole_den_old(ba, dm, 1, 0); + MultiFab e_den_old(ba, dm, 1, 0); MultiFab charge_den(ba, dm, 1, 0); MultiFab MaterialMask(ba, dm, 1, 1); MultiFab tphaseMask(ba, dm, 1, 1); @@ -160,9 +162,9 @@ void main_main (c_FerroX& rFerroX) amrex::Print() << "contains_SC = " << contains_SC << "\n"; #ifdef AMREX_USE_EB - MultiFab Plt(ba, dm, 18, 0, MFInfo(), *rGprop.pEB->p_factory_union); + MultiFab Plt(ba, dm, 24, 0, MFInfo(), *rGprop.pEB->p_factory_union); #else - MultiFab Plt(ba, dm, 18, 0); + MultiFab Plt(ba, dm, 24, 0); #endif SetPoissonBC(rFerroX, LinOpBCType_2d, all_homogeneous_boundaries, some_functionbased_inhomogeneous_boundaries, some_constant_inhomogeneous_boundaries); @@ -346,10 +348,16 @@ void main_main (c_FerroX& rFerroX) MultiFab::Copy(Plt, angle_beta, 0, 15, 1, 0); MultiFab::Copy(Plt, angle_theta, 0, 16, 1, 0); MultiFab::Copy(Plt, Phidiff, 0, 17, 1, 0); + MultiFab::Copy(Plt, Jn[0], 0, 18, 1, 0); + MultiFab::Copy(Plt, Jn[1], 0, 19, 1, 0); + MultiFab::Copy(Plt, Jn[2], 0, 20, 1, 0); + MultiFab::Copy(Plt, Jp[0], 0, 21, 1, 0); + MultiFab::Copy(Plt, Jp[1], 0, 22, 1, 0); + MultiFab::Copy(Plt, Jp[2], 0, 23, 1, 0); #ifdef AMREX_USE_EB - amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "Jnx", "Jny", "Jnz", "Jpx", "Jpy", "Jpz"}, geom, time, step); #else - amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "Jnx", "Jny", "Jnz", "Jpx", "Jpy", "Jpz"}, geom, time, step); #endif } @@ -387,6 +395,8 @@ void main_main (c_FerroX& rFerroX) err = 1.0; iter = 0; + MultiFab::Copy(e_den_old, e_den, 0, 0, 1, 0); + MultiFab::Copy(hole_den_old, hole_den, 0, 0, 1, 0); // iterate to compute Phi^{n+1,*} while(err > tol){ @@ -412,9 +422,9 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); - ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, geom, MaterialMask); + //ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -459,6 +469,8 @@ void main_main (c_FerroX& rFerroX) err = 1.0; iter = 0; + MultiFab::Copy(e_den_old, e_den, 0, 0, 1, 0); + MultiFab::Copy(hole_den_old, hole_den, 0, 0, 1, 0); // iterate to compute Phi^{n+1} while(err > tol){ @@ -484,8 +496,8 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); - ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, geom, MaterialMask); + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + //ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -568,10 +580,16 @@ void main_main (c_FerroX& rFerroX) MultiFab::Copy(Plt, angle_beta, 0, 15, 1, 0); MultiFab::Copy(Plt, angle_theta, 0, 16, 1, 0); MultiFab::Copy(Plt, Phidiff, 0, 17, 1, 0); + MultiFab::Copy(Plt, Jn[0], 0, 18, 1, 0); + MultiFab::Copy(Plt, Jn[1], 0, 19, 1, 0); + MultiFab::Copy(Plt, Jn[2], 0, 20, 1, 0); + MultiFab::Copy(Plt, Jp[0], 0, 21, 1, 0); + MultiFab::Copy(Plt, Jp[1], 0, 22, 1, 0); + MultiFab::Copy(Plt, Jp[2], 0, 23, 1, 0); #ifdef AMREX_USE_EB - amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "Jnx", "Jny", "Jnz", "Jpx", "Jpy", "Jpz"}, geom, time, step); #else - amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "Jnx", "Jny", "Jnz", "Jpx", "Jpy", "Jpz"}, geom, time, step); #endif } @@ -596,6 +614,8 @@ void main_main (c_FerroX& rFerroX) err = 1.0; iter = 0; + MultiFab::Copy(e_den_old, e_den, 0, 0, 1, 0); + MultiFab::Copy(hole_den_old, hole_den, 0, 0, 1, 0); // iterate to compute Phi^{n+1} with new Dirichlet value while(err > tol){ @@ -621,8 +641,8 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); - ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, geom, MaterialMask); + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + //ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates From 00aeeba728cf647924e11c3c730f3f251898417e Mon Sep 17 00:00:00 2001 From: prkkumar Date: Tue, 8 Aug 2023 21:32:21 -0700 Subject: [PATCH 3/5] update newton's method dF_dPhi for drift-diffusion --- Source/Solver/ElectrostaticSolver.H | 17 +++++++++++++ Source/Solver/ElectrostaticSolver.cpp | 36 +++++++++++++++++++++++++++ Source/main.cpp | 19 +++++++------- 3 files changed, 63 insertions(+), 9 deletions(-) diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 592c99e..9bd95e0 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -51,6 +51,23 @@ void dF_dPhi(MultiFab& alpha_cc, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); +void dF_dPhi_DD(MultiFab& alpha_cc, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + Array& P_old, + Array& Jn, + Array& Jp, + MultiFab& rho, + MultiFab& e_den, + MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, + MultiFab& MaterialMask, + MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi); + void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, MultiFab& alpha_cc); diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 486ca6d..8413d03 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -109,6 +109,42 @@ void dF_dPhi(MultiFab& alpha_cc, MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0); } + +void dF_dPhi_DD(MultiFab& alpha_cc, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + Array& P_old, + Array& Jn, + Array& Jp, + MultiFab& rho, + MultiFab& e_den, + MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, + MultiFab& MaterialMask, + MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi) + +{ + + MultiFab PoissonPhi_plus_delta(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 0); + MultiFab PoissonRHS_phi_plus_delta(PoissonRHS.boxArray(), PoissonRHS.DistributionMap(), 1, 0); + + MultiFab::Copy(PoissonPhi_plus_delta, PoissonPhi, 0, 0, 1, 0); + PoissonPhi_plus_delta.plus(delta, 0, 1, 0); + + // Calculate rho from Phi in SC region + //ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, rho, e_den, p_den, e_den_old, p_den_old, geom, MaterialMask); + + //Compute RHS of Poisson equation + ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); + + MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0); +} + void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, MultiFab& alpha_cc) diff --git a/Source/main.cpp b/Source/main.cpp index dacd46c..6f87c73 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -404,7 +404,8 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new_pre, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + dF_dPhi_DD(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -422,9 +423,9 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); - //ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -478,7 +479,7 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -496,8 +497,8 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); - //ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); + //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -623,7 +624,7 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -641,8 +642,8 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); - //ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); + //ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates From 221ef30e0fb7c1e7e75a0be7fe115f682dc2e7e7 Mon Sep 17 00:00:00 2001 From: prkkumar Date: Tue, 8 Aug 2023 23:02:21 -0700 Subject: [PATCH 4/5] update newtons method everywhere --- Source/main.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Source/main.cpp b/Source/main.cpp index 6f87c73..e1e40ee 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -137,13 +137,16 @@ void main_main (c_FerroX& rFerroX) } //Initialize material mask - InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); - //InitializeMaterialMask(rFerroX, geom, MaterialMask); + //InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); + InitializeMaterialMask(rFerroX, geom, MaterialMask); if(Coordinate_Transformation == 1){ Initialize_tphase_Mask(rFerroX, geom, tphaseMask); Initialize_Euler_angles(rFerroX, geom, angle_alpha, angle_beta, angle_theta); } else { tphaseMask.setVal(0.); + angle_alpha.setVal(0.); + angle_beta.setVal(0.); + angle_theta.setVal(0.); } //Solver for Poisson equation @@ -479,7 +482,8 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + dF_dPhi_DD(alpha_cc, PoissonRHS, PoissonPhi, P_new, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -624,7 +628,8 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + dF_dPhi_DD(alpha_cc, PoissonRHS, PoissonPhi, P_old, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); From 612beeb96a969f89f02e394b58492c0af64fd970 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Fri, 11 Aug 2023 11:37:18 -0700 Subject: [PATCH 5/5] add dFdPhi function with drift-diffusion --- Source/Solver/ElectrostaticSolver.H | 28 ++++++++++++++++---- Source/Solver/ElectrostaticSolver.cpp | 38 ++++++++++++++++++++++++++- Source/main.cpp | 10 +++++-- 3 files changed, 68 insertions(+), 8 deletions(-) diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 592c99e..50534db 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -39,18 +39,36 @@ void InitializePermittivity(std::array& prob_hi); void dF_dPhi(MultiFab& alpha_cc, - MultiFab& PoissonRHS, - MultiFab& PoissonPhi, - Array& P_old, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + Array& P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, + MultiFab& MaterialMask, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); +void dF_dPhi_DD(MultiFab& alpha_cc, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + Array& P_old, + Array& Jn, + Array& Jp, + MultiFab& rho, + MultiFab& e_den, + MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, + MultiFab& MaterialMask, + MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi); + + void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, MultiFab& alpha_cc); diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 486ca6d..79cedfb 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -102,13 +102,49 @@ void dF_dPhi(MultiFab& alpha_cc, PoissonPhi_plus_delta.plus(delta, 0, 1, 0); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); + ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, MaterialMask); //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0); } + +void dF_dPhi_DD(MultiFab& alpha_cc, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + Array& P_old, + Array& Jn, + Array& Jp, + MultiFab& rho, + MultiFab& e_den, + MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, + MultiFab& MaterialMask, + MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi) + +{ + + MultiFab PoissonPhi_plus_delta(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 0); + MultiFab PoissonRHS_phi_plus_delta(PoissonRHS.boxArray(), PoissonRHS.DistributionMap(), 1, 0); + + MultiFab::Copy(PoissonPhi_plus_delta, PoissonPhi, 0, 0, 1, 0); + PoissonPhi_plus_delta.plus(delta, 0, 1, 0); + + // Calculate rho from Phi in SC region + //ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi_plus_delta, Jn, Jp, rho, e_den, p_den, e_den_old, p_den_old, geom, MaterialMask); + + //Compute RHS of Poisson equation + ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); + + MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0); +} + void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, MultiFab& alpha_cc) diff --git a/Source/main.cpp b/Source/main.cpp index dacd46c..2e9f393 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -137,13 +137,16 @@ void main_main (c_FerroX& rFerroX) } //Initialize material mask - InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); - //InitializeMaterialMask(rFerroX, geom, MaterialMask); + //InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); + InitializeMaterialMask(rFerroX, geom, MaterialMask); if(Coordinate_Transformation == 1){ Initialize_tphase_Mask(rFerroX, geom, tphaseMask); Initialize_Euler_angles(rFerroX, geom, angle_alpha, angle_beta, angle_theta); } else { tphaseMask.setVal(0.); + angle_alpha.setVal(0.); + angle_beta.setVal(0.); + angle_theta.setVal(0.); } //Solver for Poisson equation @@ -405,6 +408,7 @@ void main_main (c_FerroX& rFerroX) ComputePoissonRHS(PoissonRHS, P_new_pre, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi_DD(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -479,6 +483,7 @@ void main_main (c_FerroX& rFerroX) ComputePoissonRHS(PoissonRHS, P_new, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi_DD(alpha_cc, PoissonRHS, PoissonPhi, P_new, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -624,6 +629,7 @@ void main_main (c_FerroX& rFerroX) ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi_DD(alpha_cc, PoissonRHS, PoissonPhi, P_old, Jn, Jp, charge_den, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc);