Skip to content

Commit

Permalink
first and 2nd order time integrator and nodality of gradient energy
Browse files Browse the repository at this point in the history
  • Loading branch information
@prkkumar
prkkumar committed May 15, 2024
1 parent f7131d6 commit 4744830
Show file tree
Hide file tree
Showing 3 changed files with 149 additions and 74 deletions.
109 changes: 80 additions & 29 deletions Source/FieldSolver/FerroE/Eff_Field_Grad.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,100 +6,151 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real UpwardDx (
amrex::Array4<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> 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<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> 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<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> 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<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> 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<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> 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<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> 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<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> const& fe) {
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> 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]);

}

}
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDy (
amrex::Array4<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> 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<amrex::Real, 3> dx, amrex::Array4<amrex::Real> 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];

}

}
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDz (
amrex::Array4<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> 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<amrex::Real, 3> dx, amrex::Array4<amrex::Real> 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]);
}

}
}
Expand Down
112 changes: 68 additions & 44 deletions Source/FieldSolver/FerroE/FerroE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::Real, 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 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
Expand Down Expand Up @@ -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) {
Expand All @@ -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) {
Expand All @@ -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);
}
}
);
Expand Down
2 changes: 1 addition & 1 deletion Source/Make.WarpX
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ ifeq ($(USE_LLG),TRUE)
endif

ifeq ($(USE_FERROE),TRUE)
USERSuffix := $(USERSuffix).FERRO
USERSuffix := $(USERSuffix).FERROE
DEFINES += -DWARPX_FERROE
endif

Expand Down

0 comments on commit 4744830

Please sign in to comment.