added total energy calculation

Source/Solver/Energy_Calculation.H

@@ -0,0 +1,28 @@
+#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 ComputeGradientEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
+                            MultiFab& MaterialMask,
+                           const Geometry& geom,
+                           Real g11,
+                           Real g44);
+
+Real ComputeElectrostaticEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
+                                 Array<MultiFab, AMREX_SPACEDIM>& E,
+                                const MultiFab& MaterialMask,
+                                const Geometry& geom);
+
+#endif

Source/Solver/Energy_Calculation.cpp

@@ -0,0 +1,160 @@
+#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 ComputeGradientEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
+                            MultiFab& MaterialMask,
+                           const Geometry& geom,
+                           Real g11,
+                           Real g44)
+{
+    const Real dx = geom.CellSize(0);
+    const Real dy = geom.CellSize(1);
+    const Real dz = geom.CellSize(2);
+    const Real dV = dx * dy * dz;
+
+    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);  // Assume gradient only for 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 = (pz(i+1,j,k) - pz(i-1,j,k)) / (2.0 * dx);
+                    Real dPdy = (pz(i,j+1,k) - pz(i,j-1,k)) / (2.0 * dy);
+                    Real dPdz = (pz(i,j,k+1) - pz(i,j,k-1)) / (2.0 * dz);
+
+                    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;
+}

Source/Solver/Make.package

@@ -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

Source/main.cpp

@@ -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"
@@ -258,7 +259,23 @@ 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;
+            printf("g11_1 value is:%g \n", g11);
+            printf("g44_1 value is:%g \n", g44);
+            printf("alpha value is:%g \n", alpha);
+            printf("beta value is:%g \n", beta);
+            printf("gamma value is:%g \n", FerroX::gamma);
+            printf("dx is:%g \n", geom.CellSize(0));
+            printf("Landau_Energy value is:%g \n", Landau_Energy);
+            printf("Gradient_Energy value is:%g \n", Gradient_Energy);
+            printf("Electrostatic_Energy value is:%g \n", Electrostatic_Energy);
+            printf("Total_Energy value is:%g \n", Total_Energy);
 
+
+
             // 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);
 
@@ -291,7 +308,7 @@ 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);
-
+                printf("g11_2 value is:%g", g11);
                 // 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);