diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 20ef757..acedbe7 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -184,9 +184,11 @@ int FerroX::plot_PhiDiff; // multimaterial stack geometry AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE_lo; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE1_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::FE_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE_hi; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE1_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::FE_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_hi; @@ -233,6 +235,8 @@ AMREX_GPU_MANAGED amrex::Real FerroX::lambda; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Remnant_P; +AMREX_GPU_MANAGED amrex::Real FerroX::noise_amplitude; +AMREX_GPU_MANAGED amrex::Real FerroX::theta_dep; //problem type : initialization of P for 2D/3D/convergence problems AMREX_GPU_MANAGED int FerroX::prob_type; @@ -400,6 +404,12 @@ void InitializeFerroXNamespace(const amrex::GpuArray DE_lo; + extern AMREX_GPU_MANAGED amrex::GpuArray DE1_lo; extern AMREX_GPU_MANAGED amrex::GpuArray FE_lo; extern AMREX_GPU_MANAGED amrex::GpuArray SC_lo; extern AMREX_GPU_MANAGED amrex::GpuArray DE_hi; + extern AMREX_GPU_MANAGED amrex::GpuArray DE1_hi; extern AMREX_GPU_MANAGED amrex::GpuArray FE_hi; extern AMREX_GPU_MANAGED amrex::GpuArray SC_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_hi; @@ -73,7 +75,8 @@ namespace FerroX { extern AMREX_GPU_MANAGED amrex::GpuArray P_BC_flag_lo; extern AMREX_GPU_MANAGED amrex::GpuArray P_BC_flag_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Remnant_P; - + extern AMREX_GPU_MANAGED amrex::Real noise_amplitude; + extern AMREX_GPU_MANAGED amrex::Real theta_dep; //problem type : initialization of P for 2D/3D/convergence problems extern AMREX_GPU_MANAGED int prob_type; diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H index 80ee376..4f918d0 100644 --- a/Source/Solver/DerivativeAlgorithm.H +++ b/Source/Solver/DerivativeAlgorithm.H @@ -193,8 +193,9 @@ using namespace FerroX; } else if (P_BC_flag_lo[2] == 1){ - Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); - return (dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); // dP/dz = P_lo/lambda; + //Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); + Real dF_lo = F(i,j,k)/lambda; + return (dx[2]*dF_lo - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); // dP/dz = P_lo/lambda; // Real F_lo = (9. * F(i,j,k) - F(i,j,k+1)) / (3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return -(dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]);// dP/dz = P_lo/lambda; @@ -222,8 +223,9 @@ using namespace FerroX; } else if (P_BC_flag_hi[2] == 1){ - Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); - return (dx[2]*F_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; + //Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); + Real dF_hi = -F(i,j,k)/lambda; + return (dx[2]*dF_hi + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; // Real F_hi = (9. * F(i,j,k) - F(i,j,k-1)) / ( - 3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return -(dx[2]*F_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; @@ -407,8 +409,9 @@ using namespace FerroX; } else if (P_BC_flag_lo[2] == 1){ - Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); - return (-dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];//dPdz = P_lo/lambda; + //Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); + Real dF_lo = F(i,j,k)/lambda; + return (-dx[2]*dF_lo - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];//dPdz = P_lo/lambda; // Real F_lo = (9. * F(i,j,k) - F(i,j,k+1)) / (3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return (-dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];// dPdz = P_lo/lambda; @@ -436,8 +439,9 @@ using namespace FerroX; } else if (P_BC_flag_hi[2] == 1){ - Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); - return (dx[2]*F_hi/lambda - F(i,j,k) + F(i,j,k-1))/dx[2]/dx[2];//dPdz = P_hi/lambda; + //Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); + Real dF_hi = -F(i,j,k)/lambda; + return (dx[2]*dF_hi - F(i,j,k) + F(i,j,k-1))/dx[2]/dx[2];//dPdz = P_hi/lambda; // Real F_hi = (9. * F(i,j,k) - F(i,j,k-1)) / ( - 3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return (dx[2]*F_hi/lambda - F(i,j,k) + F(i,j,k-1))/dx[2]/dx[2]; // dPdz = P_hi/lambda; diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 128e6cb..72c6b83 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -63,6 +63,7 @@ void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& alpha_cc); void SetPhiBC_z(MultiFab& PossonPhi, const amrex::GpuArray& n_cell, const Geometry& geom); +void SetPhiBC_z(MultiFab& PoissonPhi, Array &P_old, MultiFab& MaterialMask, const amrex::GpuArray& n_cell, const Geometry& geom); void SetPoissonBC(c_FerroX& rFerroX, std::array,2>& LinOpBCType_2d, bool& all_homogeneous_boundaries, bool& some_functionbased_inhomogeneous_boundaries, bool& some_constant_inhomogeneous_boundaries); @@ -70,11 +71,14 @@ void Fill_Constant_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& Poisson void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time); void CheckSteadyState(MultiFab& PoissonPhi, MultiFab& PoissonPhi_Old, MultiFab& Phidiff, Real phi_tolerance, int step, int& steady_state_step, int& inc_step); +void SetNucleation(Array &P_old, MultiFab& NucleationMask, const amrex::GpuArray& n_cell); void SetupMLMG(std::unique_ptr& pMLMG, std::unique_ptr& p_mlabec, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + Array& P_old, + MultiFab& MaterialMask, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info); #ifdef AMREX_USE_EB @@ -101,7 +105,8 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, MultiFab& MaterialMask, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& n_cell, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); #ifdef AMREX_USE_EB diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 9286cf5..9a20ece 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -3,7 +3,7 @@ #include "ChargeDensity.H" #include "Utils/eXstaticUtils/eXstaticUtil.H" #include "Utils/FerroXUtils/FerroXUtil.H" - +#include "TotalEnergyDensity.H" void ComputePoissonRHS(MultiFab& PoissonRHS, Array &P_old, @@ -12,6 +12,18 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, const Geometry& geom) { + //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); + + //Real FE_thickness = FE_hi[2] - FE_lo[2]; + //Real one_m_theta = 1.0 - theta_dep; + //Real E_dep = average_P_r/(epsilonZ_fe*epsilon_0)*one_m_theta; + + //amrex::Print() << " E_dep = " << E_dep << "\n"; + //amrex::Print() << " FE_thickness = " << FE_thickness << "\n"; for ( MFIter mfi(PoissonRHS); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); @@ -73,7 +85,7 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, } else { //mask(i,j,k) == 0.0 FE region RHS(i,j,k) = - (R_11*DPDx(pOld_p, mask, i, j, k, dx) + R_12*DPDy(pOld_p, mask, i, j, k, dx) + R_13*DPDz(pOld_p, mask, i, j, k, dx)) - (R_21*DPDx(pOld_q, mask, i, j, k, dx) + R_22*DPDy(pOld_q, mask, i, j, k, dx) + R_23*DPDz(pOld_q, mask, i, j, k, dx)) - - (R_31*DPDx(pOld_r, mask, i, j, k, dx) + R_32*DPDy(pOld_r, mask, i, j, k, dx) + R_33*DPDz(pOld_r, mask, i, j, k, dx)); + - (R_31*DPDx(pOld_r, mask, i, j, k, dx) + R_32*DPDy(pOld_r, mask, i, j, k, dx) + R_33*DPDz(pOld_r, mask, i, j, k, dx));// - E_dep*FE_thickness; } @@ -499,8 +511,47 @@ void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& Po } } +void SetPhiBC_z(MultiFab& PoissonPhi, Array &P_old, MultiFab& MaterialMask, const amrex::GpuArray& n_cell, const Geometry& geom) +{ +// +// 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); +// +// Real FE_thickness = FE_hi[2] - FE_lo[2]; +// Real one_m_theta = 1.0 - theta_dep; +// Real E_dep = -1.0*average_P_r/(epsilonZ_fe*epsilon_0)*one_m_theta; +// +// amrex::Print() << "average_P_r = " << average_P_r << ", "<< ", E_dep = " << E_dep << ", theta = " << theta_dep << ", phi_Bc_hi" << Phi_Bc_hi << "\n"; +// + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.growntilebox(1); + + const Array4& Phi = PoissonPhi.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + if(k < 0) { + Phi(i,j,k) = Phi_Bc_lo; + } else if(k >= n_cell[2]){ + amrex::Real Eg = bandgap; + amrex::Real Chi = affinity; + amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q; + amrex::Real phi_m = use_work_function ? metal_work_function : phi_ref; //in eV When not used, applied voltgae is set as the potential on the metal interface + Phi(i,j,k) = (Phi_Bc_hi - (phi_m - phi_ref))*theta_dep;// - E_dep*FE_thickness; //multiplying by theta_dep + //if(i == 5 && j == 5) amrex::Print() << "Phi(i,j,k) = "<< Phi(i,j,k) << std::endl; + } + }); + } + PoissonPhi.FillBoundary(geom.periodicity()); +} + void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray& n_cell, const Geometry& geom) { + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(1); @@ -516,13 +567,81 @@ void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray amrex::Real Chi = affinity; amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q; amrex::Real phi_m = use_work_function ? metal_work_function : phi_ref; //in eV When not used, applied voltgae is set as the potential on the metal interface - Phi(i,j,k) = Phi_Bc_hi - (phi_m - phi_ref); + Phi(i,j,k) = (Phi_Bc_hi - (phi_m - phi_ref)); + //if(i == 5 && j == 5) amrex::Print() << "Phi(i,j,k) = "<< Phi(i,j,k) << std::endl; } }); } PoissonPhi.FillBoundary(geom.periodicity()); } +void SetNucleation(Array &P_old, MultiFab& NucleationMask, const amrex::GpuArray& n_cell) +{ + int seed = random_seed; + + int nprocs = ParallelDescriptor::NProcs(); + + if (prob_type == 1) { + amrex::InitRandom(seed , nprocs, seed ); // give all MPI ranks the same seed + } else { + amrex::InitRandom(seed+ParallelDescriptor::MyProc(), nprocs, seed+ParallelDescriptor::MyProc()); // give all MPI ranks a different seed + } + + int nrand = n_cell[0]*n_cell[2]; + amrex::Gpu::ManagedVector rngs(nrand, 0.0); + + // generate random numbers on the host + for (int i=0; i &pOld_p = P_old[0].array(mfi); + const Array4 &pOld_q = P_old[1].array(mfi); + const Array4 &pOld_r = P_old[2].array(mfi); + const Array4& mask = NucleationMask.array(mfi); + + Real* rng = rngs.data(); + + amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::RandomEngine const& engine) noexcept + { + if (mask(i,j,k) == 0.) { + if (prob_type == 1) { //2D + if (rng[i + k*n_cell[2]] <= 0.02){ + pOld_p(i,j,k) = Remnant_P[0]; + pOld_q(i,j,k) = Remnant_P[1]; + pOld_r(i,j,k) = Remnant_P[2]; + } else if (rng[i + k*n_cell[2]] <= 0.04){ + pOld_p(i,j,k) = -Remnant_P[0]; + pOld_q(i,j,k) = -Remnant_P[1]; + pOld_r(i,j,k) = -Remnant_P[2]; + } else { + pOld_p(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[0]*noise_amplitude; + pOld_q(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[1]*noise_amplitude; + pOld_r(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[2]*noise_amplitude; + } + } else if (prob_type == 2) { //3D + Real rand = Random(engine); + if (rand <= 0.04) { + pOld_p(i,j,k) = Remnant_P[0]; + pOld_q(i,j,k) = Remnant_P[1]; + pOld_r(i,j,k) = Remnant_P[2]; + } else { + pOld_p(i,j,k) += (-1.0 + 2.0*rand)*Remnant_P[0]*noise_amplitude; + pOld_q(i,j,k) += (-1.0 + 2.0*rand)*Remnant_P[1]*noise_amplitude; + pOld_r(i,j,k) += (-1.0 + 2.0*rand)*Remnant_P[2]*noise_amplitude; + } + } + } + }); + } + +} + void CheckSteadyState(MultiFab& PoissonPhi, MultiFab& PoissonPhi_Old, MultiFab& Phidiff, Real phi_tolerance, int step, int& steady_state_step, int& inc_step) { @@ -564,6 +683,8 @@ void SetupMLMG(std::unique_ptr& pMLMG, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + Array& P_old, + MultiFab& MaterialMask, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info) { auto& rGprop = rFerroX.get_GeometryProperties(); @@ -597,7 +718,8 @@ void SetupMLMG(std::unique_ptr& pMLMG, PoissonPhi.FillBoundary(geom.periodicity()); // set Dirichlet BC by reading in the ghost cell values - SetPhiBC_z(PoissonPhi, n_cell, geom); + SetPhiBC_z(PoissonPhi, n_cell, geom); + //SetPhiBC_z(PoissonPhi, P_old, MaterialMask, n_cell, geom); p_mlabec->setLevelBC(amrlev, &PoissonPhi); // (A*alpha_cc - B * div beta grad) phi = rhs @@ -689,7 +811,8 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, MultiFab& MaterialMask, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& n_cell, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { @@ -714,6 +837,8 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, //Initial guess for phi PoissonPhi.setVal(0.); + SetPhiBC_z(PoissonPhi, n_cell, geom); + //SetPhiBC_z(PoissonPhi, P_old, MaterialMask, n_cell, geom); //Poisson Solve pMLMG->solve({&PoissonPhi}, {&PoissonRHS}, 1.e-10, -1); diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 447adf4..0fde6c4 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -192,6 +192,8 @@ void InitializeMaterialMask(MultiFab& MaterialMask, mask(i,j,k) = 0.; } else if (x <= DE_hi[0] && x >= DE_lo[0] && y <= DE_hi[1] && y >= DE_lo[1] && z <= DE_hi[2] && z >= DE_lo[2]) { mask(i,j,k) = 1.; + } else if (x <= DE1_hi[0] && x >= DE1_lo[0] && y <= DE1_hi[1] && y >= DE1_lo[1] && z <= DE1_hi[2] && z >= DE1_lo[2]) { + mask(i,j,k) = 1.; } else if (x <= SC_hi[0] && x >= SC_lo[0] && y <= SC_hi[1] && y >= SC_lo[1] && z <= SC_hi[2] && z >= SC_lo[2]) { mask(i,j,k) = 2.; if (x <= Channel_hi[0] && x >= Channel_lo[0] && y <= Channel_hi[1] && y >= Channel_lo[1] && z <= Channel_hi[2] && z >= Channel_lo[2]){ @@ -218,7 +220,8 @@ void InitializeMaterialMask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& M for (MFIter mfi(MaterialMask, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const auto& mask_arr = MaterialMask.array(mfi); - const auto& bx = mfi.tilebox(); + // const auto& bx = mfi.tilebox(); + const Box& bx = mfi.growntilebox(MaterialMask.nGrow()); std::string m_mask_s; std::unique_ptr m_mask_parser; diff --git a/Source/Solver/TotalEnergyDensity.H b/Source/Solver/TotalEnergyDensity.H index 15a1326..425ab3e 100644 --- a/Source/Solver/TotalEnergyDensity.H +++ b/Source/Solver/TotalEnergyDensity.H @@ -16,3 +16,9 @@ void CalculateTDGL_RHS(Array &GL_rhs, const Geometry& geom, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); + +void Compute_P_Sum(const std::array& P, Real& sum); + +void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count); + +void Compute_P_av(const std::array& 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 53ef0bc..342a297 100644 --- a/Source/Solver/TotalEnergyDensity.cpp +++ b/Source/Solver/TotalEnergyDensity.cpp @@ -14,7 +14,18 @@ void CalculateTDGL_RHS(Array &GL_rhs, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { - // loop over boxes +// 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 one_m_theta = 1.0 - theta_dep; +// Real E_dep = average_P_r/(epsilon_0*epsilonZ_fe)*one_m_theta; +// + // loop over boxes for ( MFIter mfi(P_old[0]); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); @@ -148,7 +159,7 @@ void CalculateTDGL_RHS(Array &GL_rhs, GL_RHS_r(i,j,k) = -1.0 * Gam(i,j,k) * ( dFdPr_Landau + dFdPr_grad - - Er(i,j,k) + - Er(i,j,k) //+ E_dep ); if (is_polarization_scalar == 1){ @@ -167,4 +178,73 @@ void CalculateTDGL_RHS(Array &GL_rhs, } } +void Compute_P_Sum(const std::array& P, Real& sum) +{ + + // Initialize to zero + sum = 0.; + + ReduceOps reduce_op; + + ReduceData 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 reduce_op; + + ReduceData 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& 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; +} diff --git a/Source/main.cpp b/Source/main.cpp index 930caaf..90a6322 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -145,6 +145,7 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.setVal(0.); PoissonRHS.setVal(0.); tphaseMask.setVal(0.); + MaterialMask.setVal(0.); angle_alpha.setVal(0.); angle_beta.setVal(0.); angle_theta.setVal(0.); @@ -190,7 +191,7 @@ void main_main (c_FerroX& rFerroX) int linop_maxorder = 2; int amrlev = 0; //refers to the setcoarsest level of the solve - SetupMLMG(pMLMG, p_mlabec, LinOpBCType_2d, n_cell, beta_face, rFerroX, PoissonPhi, time, info); + SetupMLMG(pMLMG, p_mlabec, LinOpBCType_2d, n_cell, beta_face, P_old, MaterialMask, rFerroX, PoissonPhi, time, info); #ifdef AMREX_USE_EB std::unique_ptr p_mlebabec; @@ -201,7 +202,8 @@ void main_main (c_FerroX& rFerroX) //InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, geom, prob_lo, prob_hi);//old InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based - + SetNucleation(P_old, MaterialMask, n_cell); + #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_old, charge_den, e_den, hole_den, MaterialMask, @@ -209,7 +211,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_old, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif // Calculate E from Phi @@ -244,6 +246,7 @@ void main_main (c_FerroX& rFerroX) P_new_pre[i].FillBoundary(geom.periodicity()); } + #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_new_pre, charge_den, e_den, hole_den, MaterialMask, @@ -251,7 +254,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_new_pre, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif if (TimeIntegratorOrder == 1) { @@ -263,6 +266,7 @@ void main_main (c_FerroX& rFerroX) P_old[i].FillBoundary(geom.periodicity()); P_new_pre[i].FillBoundary(geom.periodicity()); } + SetNucleation(P_old, MaterialMask, n_cell); } else { @@ -282,7 +286,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_new, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif // copy new solution into old solution @@ -332,6 +336,7 @@ void main_main (c_FerroX& rFerroX) amrex::Print() << "step = " << step << ", Phi_Bc_hi = " << Phi_Bc_hi << ", num_Vapp = " << num_Vapp << ", sign = " << sign << std::endl; // Set Dirichlet BC for Phi in z + //SetPhiBC_z(PoissonPhi, P_old, MaterialMask, n_cell, geom); SetPhiBC_z(PoissonPhi, n_cell, geom); // set Dirichlet BC by reading in the ghost cell values @@ -348,7 +353,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_old, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif }//end inc_step