From 47448306fd2579fc3537cb6616274bc01b07e32b Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Wed, 15 May 2024 13:19:02 -0700 Subject: [PATCH] first and 2nd order time integrator and nodality of gradient energy --- Source/FieldSolver/FerroE/Eff_Field_Grad.H | 109 ++++++++++++++------ Source/FieldSolver/FerroE/FerroE.cpp | 112 +++++++++++++-------- Source/Make.WarpX | 2 +- 3 files changed, 149 insertions(+), 74 deletions(-) diff --git a/Source/FieldSolver/FerroE/Eff_Field_Grad.H b/Source/FieldSolver/FerroE/Eff_Field_Grad.H index d1e711614..1ce93115d 100644 --- a/Source/FieldSolver/FerroE/Eff_Field_Grad.H +++ b/Source/FieldSolver/FerroE/Eff_Field_Grad.H @@ -6,60 +6,77 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real UpwardDx ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i+1,j,k) - F(i,j,k))/(dx[0]); + return (F(i+1,j,k,0) - F(i,j,k,0))/(dx[0]); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDx ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i,j,k) - F(i-1,j,k))/(dx[0]); + return (F(i,j,k,0) - F(i-1,j,k,0))/(dx[0]); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real UpwardDy ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i,j+1,k) - F(i,j,k))/(dx[1]); + return (F(i,j+1,k,0) - F(i,j,k,0))/(dx[1]); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDy ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i,j,k) - F(i,j-1,k))/(dx[1]); + return (F(i,j,k,0) - F(i,j-1,k,0))/(dx[1]); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real UpwardDz ( 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))/(dx[2]); + return (F(i,j,k+1,0) - F(i,j,k,0))/(dx[2]); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDz ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i,j,k) - F(i,j,k-1))/(dx[2]); + return (F(i,j,k,0) - F(i,j,k-1,0))/(dx[2]); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DoubleDx ( amrex::Array4 const& F, - int const i, int const j, int const k, amrex::GpuArray dx, amrex::Array4 const& fe) { + int const i, int const j, int const k, amrex::GpuArray dx, amrex::Array4 const& fe, int const ncomp = 0) { - if (fe(i,j,k) == 1 && fe(i+1,j,k) == 0) {//upper end of fe + if (ncomp == 0){ //x-component + if (fe(i,j,k) == 1 && fe(i+1,j,k) == 0) {//upper end of fe - return (0.0 - DownwardDx(F, i, j, k, dx))/(dx[0]); + return (0.0 - DownwardDx(F, i, j, k, dx))/dx[0]; - } else if (fe(i,j,k) == 1 && fe(i-1,j,k) == 0) { //lower end of fe + } else if (fe(i,j,k) == 1 && fe(i-1,j,k) == 0) { //lower end of fe - return (UpwardDx(F, i, j, k, dx) - 0.0)/(dx[0]); + return (UpwardDx(F, i, j, k, dx) - 0.0)/dx[0]; - } else {//bulk + } else {//bulk - return (UpwardDx(F, i, j, k, dx) - DownwardDx(F, i, j, k, dx))/(dx[0]); + return (UpwardDx(F, i, j, k, dx) - DownwardDx(F, i, j, k, dx))/dx[0]; + + } + } else { //y and z component + if (fe(i,j,k) == 1 && fe(i+1,j,k) == 0) {//upper end of fe + + return 0.5*(8.*F(i-1,j,k,0) - F(i-2,j,k,0) - 7.*F(i,j,k,0))/dx[0]/dx[0]; + + } else if (fe(i,j,k) == 1 && fe(i-1,j,k) == 0) { //lower end of fe + + return 0.5*(8.*F(i+1,j,k,0) - F(i+2,j,k,0) - 7.*F(i,j,k,0))/dx[0]/dx[0]; + + } else {//bulk + + return (UpwardDx(F, i, j, k, dx) - DownwardDx(F, i, j, k, dx))/(dx[0]); + + } } } @@ -67,19 +84,36 @@ static amrex::Real DoubleDx ( AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DoubleDy ( amrex::Array4 const& F, - int const i, int const j, int const k, amrex::GpuArray dx, amrex::Array4 const& fe) { - - if (fe(i,j,k) == 1 && fe(i,j+1,k) == 0) {//upper end of fe + int const i, int const j, int const k, amrex::GpuArray dx, amrex::Array4 const& fe, int const ncomp = 0) { + + if (ncomp == 1){ //y-component + if (fe(i,j,k) == 1 && fe(i,j+1,k) == 0) {//upper end of fe + + return (0.0 - DownwardDy(F, i, j, k, dx))/dx[1]; + + } else if (fe(i,j,k) == 1 && fe(i,j-1,k) == 0) { //lower end of fe + + return (UpwardDy(F, i, j, k, dx) - 0.0)/dx[1]; + + } else {//bulk - return (0.0 - DownwardDy(F, i, j, k, dx))/(dx[1]); + return (UpwardDy(F, i, j, k, dx) - DownwardDy(F, i, j, k, dx))/(dx[1]); - } else if (fe(i,j,k) == 1 && fe(i,j-1,k) == 0) { //lower end of fe + } + } else { //x and z component + if (fe(i,j,k) == 1 && fe(i,j+1,k) == 0) {//upper end of fe - return (UpwardDy(F, i, j, k, dx) - 0.0)/(dx[1]); + return 0.5*(8.*F(i,j-1,k,0) - F(i,j-2,k,0) - 7.*F(i,j,k,0))/dx[1]/dx[1]; - } else {//bulk + } else if (fe(i,j,k) == 1 && fe(i,j-1,k) == 0) { //lower end of fe - return (UpwardDy(F, i, j, k, dx) - DownwardDy(F, i, j, k, dx))/(dx[1]); + return 0.5*(8.*F(i,j+1,k,0) - F(i,j+2,k,0) - 7.*F(i,j,k,0))/dx[1]/dx[1]; + + } else {//bulk + + return (UpwardDy(F, i, j, k, dx) - DownwardDy(F, i, j, k, dx))/dx[1]; + + } } } @@ -87,19 +121,36 @@ static amrex::Real DoubleDy ( AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DoubleDz ( amrex::Array4 const& F, - int const i, int const j, int const k, amrex::GpuArray dx, amrex::Array4 const& fe) { - - if (fe(i,j,k) == 1 && fe(i,j,k+1) == 0) {//upper end of fe + int const i, int const j, int const k, amrex::GpuArray dx, amrex::Array4 const& fe, int const ncomp = 0) { + + if (ncomp == 2){ //z-component + if (fe(i,j,k) == 1 && fe(i,j,k+1) == 0) {//upper end of fe + + return (0.0 - DownwardDz(F, i, j, k, dx))/dx[2]; + + } else if (fe(i,j,k) == 1 && fe(i,j,k-1) == 0) { //lower end of fe + + return (UpwardDz(F, i, j, k, dx) - 0.0)/dx[2]; + + } else {//bulk + + return (UpwardDz(F, i, j, k, dx) - DownwardDz(F, i, j, k, dx))/dx[2]; + + } + } else { //x and y component + if (fe(i,j,k) == 1 && fe(i,j,k+1) == 0) {//upper end of fe + + return 0.5*(8.*F(i,j,k-1,0) - F(i,j,k-2,0) - 7.*F(i,j,k,0))/dx[2]/dx[2]; - return (0.0 - DownwardDz(F, i, j, k, dx))/(dx[2]); + } else if (fe(i,j,k) == 1 && fe(i,j,k-1) == 0) { //lower end of fe - } else if (fe(i,j,k) == 1 && fe(i,j,k-1) == 0) { //lower end of fe + return 0.5*(8.*F(i,j,k+1,0) - F(i,j,k+2,0) - 7.*F(i,j,k,0))/dx[2]/dx[2]; - return (UpwardDz(F, i, j, k, dx) - 0.0)/(dx[2]); + } else {//bulk - } else {//bulk + return (UpwardDz(F, i, j, k, dx) - DownwardDz(F, i, j, k, dx))/dx[2]; - return (UpwardDz(F, i, j, k, dx) - DownwardDz(F, i, j, k, dx))/(dx[2]); + } } } diff --git a/Source/FieldSolver/FerroE/FerroE.cpp b/Source/FieldSolver/FerroE/FerroE.cpp index 69911bf54..2e156d9c4 100644 --- a/Source/FieldSolver/FerroE/FerroE.cpp +++ b/Source/FieldSolver/FerroE/FerroE.cpp @@ -152,38 +152,68 @@ FerroE::InitializeFerroelectricMultiFabUsingParser ( *dv/dt = (E_eff - gamma*v)/mu *dP/dt = v */ -AMREX_GPU_DEVICE -amrex::Real func(amrex::Real E_eff, amrex::Real v, amrex::Real gamma, amrex::Real mu) -{ - return (E_eff - gamma*v) / mu; -} - +//AMREX_GPU_DEVICE +//amrex::Real func(amrex::Real E_eff, amrex::Real dp_dt, amrex::Real gamma, amrex::Real mu) +//{ +// amrex::Real d2p_dt2 = (E_eff - gamma*dp_dt) / mu; +// return d2p_dt2; +//} +// +//AMREX_GPU_DEVICE +//std::pair func(amrex::Real E_eff, amrex::Real dp_dt, amrex::Real gamma, amrex::Real mu) +//{ +// amrex::Real d2p_dt2 = (E_eff - gamma*dp_dt) / mu; +// return std::make_pair(dp_dt, d2p_dt2); +//} +// // RK4 time integrator AMREX_GPU_DEVICE -void update_v(amrex::Real& result, amrex::Real dt, amrex::Real E_eff, amrex::Real gamma, amrex::Real mu) +void update_p(amrex::Real& p, amrex::Real& dp_dt, amrex::Real dt, amrex::Real E_eff, amrex::Real gamma, amrex::Real mu) { - int use_RK4 = 1; - - amrex::Real k1, k2, k3, k4; - - // Calculate k1 - k1 = dt * func(E_eff, result, gamma, mu); - - // Calculate k2 - k2 = dt * func(E_eff, result + 0.5*k1, gamma, mu); - - // Calculate k3 - k3 = dt * func(E_eff, result + 0.5*k2, gamma, mu); - - // Calculate k4 - k4 = dt * func(E_eff, result + k3, gamma, mu); - - if (use_RK4){ - // Update result using weighted sum of ks - result += (1.0 / 6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4); - } else { - result += k1; - } + int use_forward_Euler = 0; +// +// //amrex::Real k1, k2, k3, k4; +// +// // Calculate k1 +// auto k1 = func(E_eff, dp_dt, gamma, mu); +// +// amrex::Real p_1 = p + 0.5*dt*k1.first; +// amrex::Real dp_dt_1 = dp_dt + 0.5*dt*k1.second; +// +// // Calculate k2 +// k2 = func(E_eff, dp_dt_1, gamma, mu); +// +// amrex::Real p_1 = p + 0.5*dt*k1.first; +// amrex::Real dp_dt_1 = dp_dt + 0.5*dt*k1.second; +// +// // Calculate k3 +// k3 = dt * func(E_eff, result + 0.5*k2, gamma, mu); +// +// // Calculate k4 +// k4 = dt * func(E_eff, result + k3, gamma, mu); +// +// if (use_RK4){ +// // Update result using weighted sum of ks +// result += (1.0 / 6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4); +// } else { +// result += k1; +// } + //auto k1 = func(E_eff, dp_dt, gamma, mu); + //p += dt*k1.first; + //dp_dt += dt*k1.second; + + if (use_forward_Euler) { + //Forward Euler + amrex::Real rhs = (E_eff - gamma*dp_dt) / mu; + dp_dt += dt*rhs; + p += dt*dp_dt; + } else { + //2nd Order + amrex::Real v_old = dp_dt; + amrex::Real fac = 0.5*dt*gamma/mu; + dp_dt = 1./(1. + fac) * (dt/mu*E_eff + (1. - fac)*v_old); + p += 0.5*dt*(v_old + dp_dt); + } } void @@ -227,13 +257,11 @@ FerroE::EvolveP (amrex::Real dt) } if (include_grad == 1){ - Ex_eff += G_11*DoubleDx(Px_arr, i, j, k, dx, fe_arr) + G_11*DoubleDy(Px_arr, i, j, k, dx, fe_arr) + G_11*DoubleDz(Px_arr, i, j, k, dx, fe_arr); + Ex_eff += G_11*DoubleDx(Px_arr, i, j, k, dx, fe_arr, 0) + G_11*DoubleDy(Px_arr, i, j, k, dx, fe_arr, 0) + G_11*DoubleDz(Px_arr, i, j, k, dx, fe_arr, 0); } - //get dPx/dt using numerical integration - update_v(Px_arr(i,j,k,1), dt, Ex_eff, gamma, mu); - //get Px - Px_arr(i,j,k,0) += dt*Px_arr(i,j,k,1); + //get Px and dPx/dt using numerical integration + update_p(Px_arr(i,j,k,0), Px_arr(i,j,k,1), dt, Ex_eff, gamma, mu); } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -246,13 +274,11 @@ FerroE::EvolveP (amrex::Real dt) } if (include_grad == 1){ - Ey_eff += G_11*DoubleDx(Py_arr, i, j, k, dx, fe_arr) + G_11*DoubleDy(Py_arr, i, j, k, dx, fe_arr) + G_11*DoubleDz(Py_arr, i, j, k, dx, fe_arr); + Ey_eff += G_11*DoubleDx(Py_arr, i, j, k, dx, fe_arr, 1) + G_11*DoubleDy(Py_arr, i, j, k, dx, fe_arr, 1) + G_11*DoubleDz(Py_arr, i, j, k, dx, fe_arr, 1); } - //get dPy/dt using numerical integration - update_v(Py_arr(i,j,k,1), dt, Ey_eff, gamma, mu); - //get Py - Py_arr(i,j,k,0) += dt*Py_arr(i,j,k,1); + //get Py and dPy/dt using numerical integration + update_p(Py_arr(i,j,k,0), Py_arr(i,j,k,1), dt, Ey_eff, gamma, mu); } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -265,13 +291,11 @@ FerroE::EvolveP (amrex::Real dt) } if (include_grad == 1){ - Ez_eff += G_11*DoubleDx(Pz_arr, i, j, k, dx, fe_arr) + G_11*DoubleDy(Pz_arr, i, j, k, dx, fe_arr) + G_11*DoubleDz(Pz_arr, i, j, k, dx, fe_arr); + Ez_eff += G_11*DoubleDx(Pz_arr, i, j, k, dx, fe_arr, 2) + G_11*DoubleDy(Pz_arr, i, j, k, dx, fe_arr, 2) + G_11*DoubleDz(Pz_arr, i, j, k, dx, fe_arr, 2); } - //get dPz/dt using numerical integration - update_v(Pz_arr(i,j,k,1), dt, Ez_eff, gamma, mu); - //get Pz - Pz_arr(i,j,k,0) += dt*Pz_arr(i,j,k,1); + //get Pz and dPz/dt using numerical integration + update_p(Pz_arr(i,j,k,0), Pz_arr(i,j,k,1), dt, Ez_eff, gamma, mu); } } ); diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 31d405d6f..c4ac76bbc 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -74,7 +74,7 @@ ifeq ($(USE_LLG),TRUE) endif ifeq ($(USE_FERROE),TRUE) - USERSuffix := $(USERSuffix).FERRO + USERSuffix := $(USERSuffix).FERROE DEFINES += -DWARPX_FERROE endif