Skip to content

added total energy calculation #67

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: development
Choose a base branch
from
Open

added total energy calculation #67

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a6a9d48
added total energy calculation
tyh0123 May 26, 2025
9222f10
fiexd gradient energy calculation issue
tyh0123 Jun 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 33 additions & 33 deletions Source/Solver/DerivativeAlgorithm.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@ using namespace FerroX;
* Perform first derivative dphi/dz */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DphiDz (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& F,
amrex::Real const z_hi, amrex::Real const z_lo,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>const& prob_lo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>const& prob_hi) {

Expand All @@ -25,27 +25,27 @@ using namespace FerroX;
* Perform first derivative dF/dx */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DFDx (
amrex::Array4<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
amrex::Array4<const amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx) {
return (F(i+1,j,k) - F(i-1,j,k))/(2.*dx[0]);
}

/**
* Perform first derivative dF/dy */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DFDy (
amrex::Array4<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
amrex::Array4<const amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx) {
return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]);
}

/**
* Perform first derivative dP/dx */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DPDx (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
) {
if (mask(i-1,j,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary

Expand Down Expand Up @@ -111,9 +111,9 @@ using namespace FerroX;
* Perform first derivative dP/dy */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DPDy (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
) {

if (mask(i,j-1,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
Expand Down Expand Up @@ -180,9 +180,9 @@ using namespace FerroX;
* Perform first derivative dP/dz */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DPDz (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
) {

if (mask(i,j,k-1) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
Expand Down Expand Up @@ -255,9 +255,9 @@ using namespace FerroX;
* Perform double derivative (d^2)P/dx^2 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDPDx (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx
) {

if (mask(i-1,j,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
Expand Down Expand Up @@ -325,9 +325,9 @@ using namespace FerroX;
* Perform double derivative (d^2)P/dy^2 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDPDy (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
) {

if (mask(i,j-1,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
Expand Down Expand Up @@ -394,9 +394,9 @@ using namespace FerroX;
* Perform double derivative (d^2)P/dz^2 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDPDz (
amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
) {

if (mask(i,j,k-1) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
Expand Down Expand Up @@ -469,29 +469,29 @@ using namespace FerroX;
/**
* Perform double derivative (d^2)P/dxdy */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDPDxDy (amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx)
static amrex::Real DoubleDPDxDy (amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx)
{
return (DPDy(F, mask, i+1, j, k, dx) - DPDy(F, mask, i-1, j, k, dx)) / 2. /dx[0];
}

/**
* Perform double derivative (d^2)P/dxdz */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDPDxDz (amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx)
static amrex::Real DoubleDPDxDz (amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx)
{
return (DPDz(F, mask, i+1, j, k, dx) - DPDz(F, mask, i-1, j, k, dx)) / 2. /dx[0];
}

/**
* Perform double derivative (d^2)P/dydz */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DoubleDPDyDz (amrex::Array4<amrex::Real> const& F,
amrex::Array4<amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx)
static amrex::Real DoubleDPDyDz (amrex::Array4<const amrex::Real> const& F,
amrex::Array4<const amrex::Real> const& mask,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx)
{
return (DPDz(F, mask, i, j+1, k, dx) - DPDz(F, mask, i, j-1, k, dx)) / 2. /dx[1];
}
Expand Down
29 changes: 29 additions & 0 deletions Source/Solver/Energy_Calculation.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
#ifndef FERROX_TOTALFREEENERGY_H_
#define FERROX_TOTALFREEENERGY_H_

#include <AMReX_Array.H>
#include <AMReX_MultiFab.H>
#include <AMReX_Geometry.H>
#include <AMReX_REAL.H>

// Computes the Landau free energy over the ferroelectric region (MaterialMask == 0.0)
Real ComputeLandauEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
MultiFab& MaterialMask,
const Geometry& geom,
Real alpha,
Real beta,
Real gamma_coeff);


Real ComputeElectrostaticEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
Array<MultiFab, AMREX_SPACEDIM>& E,
const MultiFab& MaterialMask,
const Geometry& geom);


Real ComputeGradientEnergy(const Array<MultiFab, AMREX_SPACEDIM>& P,
const MultiFab& MaterialMask,
const Geometry& geom,
Real g11,
Real g44);
#endif
121 changes: 121 additions & 0 deletions Source/Solver/Energy_Calculation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@
#include <AMReX_MultiFab.H>
#include <AMReX_Geometry.H>
#include <AMReX_REAL.H>
#include <cmath>
#include "DerivativeAlgorithm.H"
#include "Energy_Calculation.H"

using namespace amrex;

Real ComputeLandauEnergy(Array<MultiFab, AMREX_SPACEDIM>& P,
MultiFab& MaterialMask,
const Geometry& geom,
Real alpha,
Real beta,
Real gamma_coeff)
{
const Real dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);

ReduceOps<ReduceOpSum> reduce_op;
ReduceData<Real> reduce_data(reduce_op);
using ReduceTuple = typename decltype(reduce_data)::Type;

for (MFIter mfi(P[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
const Box& bx = mfi.tilebox();
auto const& px = P[0].array(mfi);
auto const& py = P[1].array(mfi);
auto const& pz = P[2].array(mfi);
auto const& mask = MaterialMask.array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE(int i, int j, int k) -> ReduceTuple {
if (mask(i,j,k) == 0.0) {
Real P2 = px(i,j,k)*px(i,j,k) + py(i,j,k)*py(i,j,k) + pz(i,j,k)*pz(i,j,k);
Real P4 = P2 * P2;
Real P6 = P4 * P2;
Real f_landau = 0.5 * alpha * P2 + 0.25 * beta * P4 + (1.0/6.0) * gamma_coeff * P6;
return {f_landau * dV};
} else {
return {0.0};
}
});
}

Real landau_energy = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(landau_energy);
return landau_energy;
}


Real ComputeGradientEnergy(const Array<MultiFab, AMREX_SPACEDIM>& P,
const MultiFab& MaterialMask,
const Geometry& geom,
Real g11,
Real g44)
{
const Real dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
const auto dx = geom.CellSizeArray(); // GpuArray<Real, 3>

ReduceOps<ReduceOpSum> reduce_op;
ReduceData<Real> 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();
auto const& pz = P[2].array(mfi); // Scalar polarization (Pz)
auto const& mask = MaterialMask.array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE(int i, int j, int k) -> ReduceTuple {
if (mask(i,j,k) == 0.0) {
Real dPdx = DPDx(pz, mask, i, j, k, dx);
Real dPdy = DPDy(pz, mask, i, j, k, dx);
Real dPdz = DPDz(pz, mask, i, j, k, dx);

Real f_grad = 0.5 * (g44 * (dPdx*dPdx + dPdy*dPdy) + g11 * dPdz*dPdz);
return {f_grad * dV};
} else {
return {0.0};
}
});
}

Real grad_energy = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(grad_energy);
return grad_energy;
}



Real ComputeElectrostaticEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
Array<MultiFab, AMREX_SPACEDIM>& E,
const MultiFab& MaterialMask,
const Geometry& geom)
{
const Real dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);

ReduceOps<ReduceOpSum> reduce_op;
ReduceData<Real> 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();
auto const& pz = P[2].array(mfi); // scalar polarization: P = Pz
auto const& ez = E[2].array(mfi); // scalar electric field: E = Ez
auto const& mask = MaterialMask.array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE(int i, int j, int k) -> ReduceTuple {
if (mask(i,j,k) == 0.0) {
Real f_elec = ez(i,j,k) * pz(i,j,k); // E · P
return {f_elec * dV};
} else {
return {0.0};
}
});
}

Real electrostatic_energy = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(electrostatic_energy);
return electrostatic_energy;
}
2 changes: 2 additions & 0 deletions Source/Solver/Make.package
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@ CEXE_sources += ElectrostaticSolver.cpp
CEXE_sources += Initialization.cpp
CEXE_sources += ChargeDensity.cpp
CEXE_sources += TotalEnergyDensity.cpp
CEXE_sources += Energy_Calculation.cpp

CEXE_headers += ElectrostaticSolver.H
CEXE_headers += Initialization.H
CEXE_headers += ChargeDensity.H
CEXE_headers += TotalEnergyDensity.H
CEXE_headers += Energy_Calculation.H

VPATH_LOCATIONS += $(CODE_HOME)/Source/Solver
INCLUDE_LOCATIONS += $(CODE_HOME)/Source/Solver
Expand Down
13 changes: 12 additions & 1 deletion Source/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "Solver/Initialization.H"
#include "Solver/ChargeDensity.H"
#include "Solver/TotalEnergyDensity.H"
#include "Solver/Energy_Calculation.H"
#include "Input/BoundaryConditions/BoundaryConditions.H"
#include "Input/GeometryProperties/GeometryProperties.H"
#include "Utils/SelectWarpXUtils/WarpXUtil.H"
Expand Down Expand Up @@ -258,7 +259,18 @@ void main_main (c_FerroX& rFerroX)

// Calculate E from Phi
ComputeEfromPhi(PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);
Real Landau_Energy = ComputeLandauEnergy(P_old, MaterialMask, geom, alpha, beta, FerroX::gamma);
Real Gradient_Energy = ComputeGradientEnergy(P_old, MaterialMask, geom, g11, g44);
Real Electrostatic_Energy = ComputeElectrostaticEnergy(P_old, E, MaterialMask, geom);
Real Total_Energy = Landau_Energy + Gradient_Energy + Electrostatic_Energy;

amrex::Print() << "Landau_Energy value is: " << Landau_Energy << '\n';
amrex::Print() << "Gradient_Energy value is: " << Gradient_Energy << '\n';
amrex::Print() << "Electrostatic_Energy value is: " << Electrostatic_Energy << '\n';
amrex::Print() << "Total_Energy value is: " << Total_Energy << '\n';



// compute f^n = f(P^n,Phi^n)
CalculateTDGL_RHS(GL_rhs, P_old, E, Gamma, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);

Expand Down Expand Up @@ -291,7 +303,6 @@ void main_main (c_FerroX& rFerroX)

//update E using PoissonPhi computed with P_new_pre
ComputeEfromPhi(PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);

// compute f^{n+1,*} = f(P^{n+1,*},Phi^{n+1,*})
CalculateTDGL_RHS(GL_rhs_pre, P_new_pre, E, Gamma, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);

Expand Down