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, + 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 f418a6f..dc9a261 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,177 @@ 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, + MultiFab& e_den_old, + MultiFab& p_den_old, + 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& 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); + 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) = 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) { + 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/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 a9e6740..e1e40ee 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); @@ -125,14 +127,26 @@ 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); + //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 @@ -151,9 +165,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); @@ -337,10 +351,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 } @@ -378,6 +398,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){ @@ -385,7 +407,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); @@ -403,7 +426,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); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -448,6 +473,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){ @@ -455,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, 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); + 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); @@ -473,7 +501,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, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -556,10 +585,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 } @@ -584,6 +619,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){ @@ -591,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, 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); + 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); @@ -609,7 +647,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, e_den_old, hole_den_old, geom, MaterialMask); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates