From 1222457cad5d784b4ac58082d72f1374484584ad Mon Sep 17 00:00:00 2001 From: Matthew Mehrtens <12023414+mcmehrtens@users.noreply.github.com> Date: Wed, 30 Sep 2026 13:17:40 -0500 Subject: [PATCH 1/4] Add singly degenerate Cahn-Hilliard mobility --- src/Integrator/CahnHilliard.H | 30 ++- src/Integrator/CahnHilliard.cpp | 178 ++++++++++++- src/Operator/Spectral/FFT.H | 26 ++ tests/SinglyDegenerateCahnHilliard/Readme.rst | 22 ++ tests/SinglyDegenerateCahnHilliard/input | 236 ++++++++++++++++++ tests/SinglyDegenerateCahnHilliard/test | 84 +++++++ 6 files changed, 572 insertions(+), 4 deletions(-) create mode 100644 tests/SinglyDegenerateCahnHilliard/Readme.rst create mode 100644 tests/SinglyDegenerateCahnHilliard/input create mode 100755 tests/SinglyDegenerateCahnHilliard/test diff --git a/src/Integrator/CahnHilliard.H b/src/Integrator/CahnHilliard.H index 0bdc85b9f7..ad2a6897d0 100644 --- a/src/Integrator/CahnHilliard.H +++ b/src/Integrator/CahnHilliard.H @@ -14,9 +14,26 @@ // // \frac{\partial\eta}{\partial t} = L\nabla^2\Big(\eta^3 - \eta - \gamma\nabla^2\eta\Big) // -// which is integrated using a forward Euler scheme. +// Constant mobility is integrated using forward Euler in realspace or a +// semi-implicit spectral scheme. With :code:`mobility = singly_degenerate`, +// the equation is instead // -// This is tested in :ref:`CahnHilliard` +// .. math:: +// +// \partial_t\eta = \nabla\cdot(M(\eta)\nabla\mu),\qquad +// M(\eta) = L\alpha^2(1-\alpha)^2 + M_0,\qquad +// \alpha = \mathrm{clip}((\eta+1)/2,0,1), +// +// where :math:`\mu=\eta^3-\eta-\gamma\nabla^2\eta` and +// :math:`M_0` is :code:`mobility_floor` (default zero). +// The realspace method uses conservative face fluxes; the spectral method +// stabilizes the flux divergence with :code:`spectral_stabilization`, default +// :math:`\gamma(L/16+M_0)`. SDCH requires periodic boundaries. Its spectral +// method assembles a uniform finest-level grid and restricts the update to +// the AMR hierarchy; realspace AMR has no flux refluxing at coarse/fine interfaces. +// The solution is not clipped, so timesteps must resolve the interface dynamics. +// +// This is tested in :ref:`CahnHilliard` and :ref:`SinglyDegenerateCahnHilliard` #ifndef INTEGRATOR_CAHNHILLIARD_H #define INTEGRATOR_CAHNHILLIARD_H @@ -60,6 +77,10 @@ protected: void AdvanceSpectral (int lev, Set::Scalar time, Set::Scalar dt); /// realspace version of advance void AdvanceReal (int lev, Set::Scalar time, Set::Scalar dt); + /// Conservative realspace update with variable mobility + void AdvanceDegenerateReal (int lev, Set::Scalar time, Set::Scalar dt); + /// Stabilized spectral update with variable mobility + void AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar dt); /// Mark any cells that need to be refined void TagCellsForRefinement (int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow) override; @@ -76,7 +97,10 @@ private: Set::Scalar L = NAN; Set::Scalar refinement_threshold = NAN; - std::string method; + std::string method; + std::string mobility = "constant"; + Set::Scalar mobility_floor = 0.0; + Set::Scalar spectral_stabilization = NAN; }; } #endif diff --git a/src/Integrator/CahnHilliard.cpp b/src/Integrator/CahnHilliard.cpp index 2a8e516ef6..abe31b6dae 100644 --- a/src/Integrator/CahnHilliard.cpp +++ b/src/Integrator/CahnHilliard.cpp @@ -1,4 +1,8 @@ #include +#include +#include + +#include "IC/Expression.H" #include "CahnHilliard.H" #include "BC/Constant.H" @@ -25,17 +29,35 @@ void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp) pp.query_default("gamma",value.gamma, 0.0005); // Mobility pp.query_default("L", value.L, 1.0); + // Mobility model (constant preserves the original update). + pp.query_validate("mobility", value.mobility, {"constant", "singly_degenerate"}); + // Additive mobility in the pure phases. + pp.query_default("mobility_floor", value.mobility_floor, 0.0); + // Coefficient of the spectral fourth-order stabilizer. + pp.query_default("spectral_stabilization", value.spectral_stabilization, + value.gamma * (value.L / 16.0 + value.mobility_floor)); // Regridding criterion pp.query_default("refinement_threshold",value.refinement_threshold, 1E100); // initial condition for :math:`\eta` - pp.select_default("eta.ic", value.ic, pp.forward_args(value.geom)); + pp.select_default("eta.ic", value.ic, pp.forward_args(value.geom)); // boundary condition for :math:`\eta` pp.select_default("eta.bc", value.bc, pp.forward_args(1)); // Which method to use - realspace or spectral method. pp.query_validate("method",value.method,{"realspace","spectral"}); + if (value.mobility == "singly_degenerate" && !pp.InTraversalMode()) + { + if (!std::isfinite(value.L) || value.L < 0.0 || + !std::isfinite(value.gamma) || value.gamma < 0.0 || + !std::isfinite(value.mobility_floor) || value.mobility_floor < 0.0 || + !std::isfinite(value.spectral_stabilization) || value.spectral_stabilization < 0.0) + Util::ParmParseException(INFO, "mobility", "SDCH requires finite nonnegative L, gamma, mobility_floor and spectral_stabilization"); + if (!value.geom[0].isAllPeriodic()) + Util::ParmParseException(INFO, "mobility", "SDCH requires periodic boundaries"); + } + value.RegisterNewFab(value.etanew_mf, value.bc, 1, 1, "eta",true); value.RegisterNewFab(value.intermediate, value.bc, 1, 1, "int",true); @@ -46,6 +68,14 @@ void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp) void CahnHilliard::Advance(int lev, Set::Scalar time, Set::Scalar dt) { + if (mobility == "singly_degenerate") + { + if (method == "realspace") + AdvanceDegenerateReal(lev, time, dt); + else if (lev == finest_level) + AdvanceDegenerateSpectral(lev, time, dt); + return; + } if (method == "realspace") AdvanceReal(lev, time, dt); else if (method == "spectral") @@ -91,6 +121,152 @@ CahnHilliard::AdvanceReal (int lev, Set::Scalar /*time*/, Set::Scalar dt) } } +namespace +{ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +Set::Scalar DegenerateMobility(Set::Scalar eta, Set::Scalar L, Set::Scalar floor) +{ + Set::Scalar alpha = std::max(0.0, std::min(1.0, 0.5 * (eta + 1.0))); + return L * alpha * alpha * (1.0 - alpha) * (1.0 - alpha) + floor; +} + +void ChemicalPotential(const amrex::MultiFab& eta_mf, amrex::MultiFab& mu_mf, + const amrex::Geometry& geometry, Set::Scalar gamma) +{ + const auto dx = geometry.CellSizeArray(); + for (amrex::MFIter mfi(eta_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto eta = eta_mf.const_array(mfi); + const auto mu = mu_mf.array(mfi); + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + mu(i,j,k) = eta(i,j,k)*eta(i,j,k)*eta(i,j,k) - eta(i,j,k) + - gamma * Numeric::Laplacian(eta, i, j, k, 0, dx.data()); + }); + } +} +} + +void +CahnHilliard::AdvanceDegenerateReal (int lev, Set::Scalar time, Set::Scalar dt) +{ + std::swap(etaold_mf[lev], etanew_mf[lev]); + etaold_mf[lev]->FillBoundary(geom[lev].periodicity()); + ChemicalPotential(*etaold_mf[lev], *intermediate[lev], geom[lev], gamma); + intermediate[lev]->FillBoundary(geom[lev].periodicity()); + if (lev > 0) + { + amrex::Vector coarse{intermediate[lev-1].get()}; + amrex::Vector fine{intermediate[lev].get()}; + amrex::Vector times{time}; + bc->define(geom[lev]); + amrex::Vector bcs{bc->GetBCRec()}; + amrex::FillPatchTwoLevels(*intermediate[lev], time, coarse, times, fine, times, + 0, 0, 1, geom[lev-1], geom[lev], *bc, 0, *bc, 0, refRatio(lev-1), + &amrex::cell_cons_interp, bcs, 0); + } + + const auto dx = geom[lev].CellSizeArray(); + const Set::Scalar L_local = L; + const Set::Scalar floor = mobility_floor; + for (amrex::MFIter mfi(*etanew_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto eta = etaold_mf[lev]->const_array(mfi); + const auto mu = intermediate[lev]->const_array(mfi); + const auto next = etanew_mf[lev]->array(mfi); + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + const Set::Scalar M = DegenerateMobility(eta(i,j,k), L_local, floor); + Set::Scalar rhs = 0.0; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) + { + const int di = dir == 0, dj = dir == 1, dk = dir == 2; + const Set::Scalar plus = 0.5 * (M + DegenerateMobility(eta(i+di,j+dj,k+dk), L_local, floor)); + const Set::Scalar minus = 0.5 * (M + DegenerateMobility(eta(i-di,j-dj,k-dk), L_local, floor)); + rhs += (plus * (mu(i+di,j+dj,k+dk) - mu(i,j,k)) + - minus * (mu(i,j,k) - mu(i-di,j-dj,k-dk))) / (dx[dir] * dx[dir]); + } + next(i,j,k) = eta(i,j,k) + dt * rhs; + }); + } + etanew_mf[lev]->FillBoundary(geom[lev].periodicity()); +} + +#ifdef ALAMO_FFT +void +CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar dt) +{ + using Operator::Spectral::FFT; + FFT fft(geom[lev]); + amrex::Vector hierarchy(lev + 1); + for (int ilev = 0; ilev <= lev; ++ilev) hierarchy[ilev] = etanew_mf[ilev].get(); + // Form derivatives on a complete grid so coarse/fine patch edges have valid data. + auto eta_mf = FFT::CompositeToUniform(hierarchy, geom, refRatio(), lev, time, *bc, 1, 1); + eta_mf->FillBoundary(geom[lev].periodicity()); + amrex::MultiFab mu_mf(eta_mf->boxArray(), eta_mf->DistributionMap(), 1, 1); + ChemicalPotential(*eta_mf, mu_mf, geom[lev], gamma); + mu_mf.FillBoundary(geom[lev].periodicity()); + amrex::Vector potential_hierarchy(lev + 1); + for (int ilev = 0; ilev <= lev; ++ilev) potential_hierarchy[ilev] = intermediate[ilev].get(); + FFT::UniformToComposite(mu_mf, potential_hierarchy, geom, refRatio(), lev, 1); + amrex::MultiFab flux_mf(eta_mf->boxArray(), eta_mf->DistributionMap(), AMREX_SPACEDIM, 0); + const auto dx = geom[lev].CellSizeArray(); + const Set::Scalar L_local = L; + const Set::Scalar floor = mobility_floor; + for (amrex::MFIter mfi(*eta_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto eta = eta_mf->const_array(mfi); + const auto mu = mu_mf.const_array(mfi); + const auto flux = flux_mf.array(mfi); + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + Set::Scalar M = DegenerateMobility(eta(i,j,k), L_local, floor); + Set::Vector grad_mu = Numeric::Gradient(mu, i, j, k, 0, dx.data()); + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) flux(i,j,k,dir) = M * grad_mu(dir); + }); + } + auto eta_hat_mf = fft.MakeSpectralFab(); + auto flux_hat_mf = fft.MakeSpectralFab(AMREX_SPACEDIM); + fft.Forward(*eta_mf, eta_hat_mf); + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) fft.Forward(flux_mf, flux_hat_mf, dir, dir); + const Set::Scalar S = spectral_stabilization; + const Set::Complex I(0.0, 1.0); + const auto length = geom[lev].Domain().length3d(); + for (amrex::MFIter mfi(eta_hat_mf, false); mfi.isValid(); ++mfi) + { + const auto eta_hat = eta_hat_mf.array(mfi); + const auto flux_hat = flux_hat_mf.const_array(mfi); + fft.ParallelForModes(mfi.tilebox(), [=] AMREX_GPU_DEVICE( + int m, int n, int p, Set::Scalar k1, Set::Scalar k2, Set::Scalar k3, Set::Scalar omega2) + { + amrex::ignore_unused(k2, k3); + // A real Fourier derivative vanishes at even-grid Nyquist modes. + if (length[0] % 2 == 0 && m == length[0]/2) k1 = 0.0; +#if AMREX_SPACEDIM >= 2 + if (length[1] % 2 == 0 && n == length[1]/2) k2 = 0.0; +#endif +#if AMREX_SPACEDIM >= 3 + if (length[2] % 2 == 0 && p == length[2]/2) k3 = 0.0; +#endif + Set::Complex div_flux = AMREX_D_TERM(I * k1 * flux_hat(m,n,p,0), + + I * k2 * flux_hat(m,n,p,1), + + I * k3 * flux_hat(m,n,p,2)); + // Keeping the zero mode unchanged conserves the mean without clipping. + eta_hat(m,n,p) += dt * div_flux / (1.0 + dt * S * omega2 * omega2); + }); + } + fft.Backward(eta_hat_mf, *eta_mf); + FFT::UniformToComposite(*eta_mf, hierarchy, geom, refRatio(), lev, 1); + for (int ilev = 0; ilev <= lev; ++ilev) etanew_mf[ilev]->FillBoundary(geom[ilev].periodicity()); +} +#else +void +CahnHilliard::AdvanceDegenerateSpectral (int, Set::Scalar, Set::Scalar) +{ + Util::Abort(INFO,"Alamo must be compiled with fft"); +} +#endif + #ifdef ALAMO_FFT void CahnHilliard::AdvanceSpectral (int lev, Set::Scalar time, Set::Scalar dt) diff --git a/src/Operator/Spectral/FFT.H b/src/Operator/Spectral/FFT.H index c7210f70af..dafa0546d9 100644 --- a/src/Operator/Spectral/FFT.H +++ b/src/Operator/Spectral/FFT.H @@ -119,6 +119,32 @@ public: }); } + template + void ParallelForModes (const amrex::Box& bx, F&& f) const + { + int length[AMREX_SPACEDIM]; + Set::Scalar kscale[AMREX_SPACEDIM]; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) + { + length[dir] = m_length[dir]; + kscale[dir] = m_kscale[dir]; + } + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int m, int n, int p) { + Set::Scalar k1 = K(0, m, length, kscale); + Set::Scalar k2 = 0.0; + Set::Scalar k3 = 0.0; +#if AMREX_SPACEDIM >= 2 + k2 = (n <= (length[1]-1)/2 ? n : n-length[1]) * kscale[1]; +#endif +#if AMREX_SPACEDIM >= 3 + k3 = (p <= (length[2]-1)/2 ? p : p-length[2]) * kscale[2]; +#endif + Set::Scalar omega2 = k1*k1 + k2*k2 + k3*k3; + f(m, n, p, k1, k2, k3, omega2); + }); + } + void Forward (const amrex::MultiFab& in, amrex::FabArray >& out, int scomp = 0, int dcomp = 0) { m_fft.forward(in, out, scomp, dcomp); diff --git a/tests/SinglyDegenerateCahnHilliard/Readme.rst b/tests/SinglyDegenerateCahnHilliard/Readme.rst new file mode 100644 index 0000000000..645f1dbbbe --- /dev/null +++ b/tests/SinglyDegenerateCahnHilliard/Readme.rst @@ -0,0 +1,22 @@ +This problem tests singly degenerate mobility in :ref:`Integrator::CahnHilliard`. + +Use :code:`mobility = singly_degenerate` with :code:`method = realspace` or +:code:`method = spectral` (requires a build configured with :code:`--fft`). +The default :code:`mobility_floor = 0` makes the pure phases immobile; +:code:`spectral_stabilization` defaults to :math:`\gamma(L/16+M_0)`. +Constant mobility remains the default. Periodic boundary conditions are required. The realspace timestep must satisfy +the explicit fourth-order stability constraint; neither method clips the solution. + +The regression sections check a linearized Fourier mode on anisotropic domains, +nonzero mobility floors, stationary phases, conservation, and energy dissipation +in two and three dimensions. MPI sections exercise patch/rank boundaries, and +the spectral AMR section exercises a partially refined grid with regridding. +Realspace AMR does not reflux coarse/fine fluxes; use a uniform grid when exact +mass conservation is required. Spectral AMR assembles a uniform finest-level +grid, which requires memory proportional to that entire grid. + +For example, run the two-dimensional serial sections with: + +.. code-block:: bash + + scripts/runtests.py --dim=2 --serial --fft tests/SinglyDegenerateCahnHilliard diff --git a/tests/SinglyDegenerateCahnHilliard/input b/tests/SinglyDegenerateCahnHilliard/input new file mode 100644 index 0000000000..6f3a3a95b5 --- /dev/null +++ b/tests/SinglyDegenerateCahnHilliard/input @@ -0,0 +1,236 @@ +#@ [2d-realspace-amr] +#@ dim = 2 +#@ check-file = amr-real +#@ args = method=realspace timestep=0.000001 amr.max_level=1 amr.nsubsteps=1 +#@ args = refinement_threshold=0.015 amr.regrid_int=2 max_step=8 +#@ args = eta.ic.expression.region0="0.15+0.1*exp(-40*((x-0.5)^2+(y-0.75)^2))" +#@ +#@ [2d-spectral-linear-odd] +#@ dim = 2 +#@ fft = true +#@ check-file = linear-odd +#@ args = method=spectral amr.n_cell=31 31 31 amr.blocking_factor=1 max_step=1 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(30*pi*y/1.5)" +#@ +#@ [2d-spectral-unstabilized] +#@ dim = 2 +#@ fft = true +#@ check-file = nonlinear +#@ args = method=spectral spectral_stabilization=0 max_step=40 +#@ +#@ [2d-realspace-linear] +#@ dim = 2 +#@ check-file = linear +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" +#@ args = max_step=1 +#@ +#@ [2d-realspace-linear-floor] +#@ dim = 2 +#@ check-file = linear-floor +#@ args = method=realspace +#@ args = mobility_floor=0.0125 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" +#@ args = max_step=1 +#@ +#@ [2d-realspace-nonlinear] +#@ dim = 2 +#@ check-file = nonlinear +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = max_step=40 +#@ +#@ [2d-realspace-stationary] +#@ dim = 2 +#@ check-file = stationary +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="-1" +#@ args = max_step=4 +#@ +#@ [2d-realspace-stationary-outside] +#@ dim = 2 +#@ check-file = stationary-outside +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="1.2" +#@ args = max_step=4 +#@ +#@ [2d-spectral-linear] +#@ dim = 2 +#@ check-file = linear +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" +#@ args = max_step=1 +#@ +#@ [2d-spectral-linear-floor] +#@ dim = 2 +#@ check-file = linear-floor +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0125 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" +#@ args = max_step=1 +#@ +#@ [2d-spectral-nonlinear] +#@ dim = 2 +#@ check-file = nonlinear +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = max_step=40 +#@ +#@ [2d-spectral-stationary] +#@ dim = 2 +#@ check-file = stationary +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="-1" +#@ args = max_step=4 +#@ +#@ [2d-spectral-stationary-outside] +#@ dim = 2 +#@ check-file = stationary-outside +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="1.2" +#@ args = max_step=4 +#@ +#@ [3d-realspace-linear] +#@ dim = 3 +#@ check-file = linear +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" +#@ args = max_step=1 +#@ +#@ [3d-realspace-linear-floor] +#@ dim = 3 +#@ check-file = linear-floor +#@ args = method=realspace +#@ args = mobility_floor=0.0125 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" +#@ args = max_step=1 +#@ +#@ [3d-realspace-nonlinear] +#@ dim = 3 +#@ check-file = nonlinear +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = max_step=40 +#@ +#@ [3d-realspace-stationary] +#@ dim = 3 +#@ check-file = stationary +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="-1" +#@ args = max_step=4 +#@ +#@ [3d-realspace-stationary-outside] +#@ dim = 3 +#@ check-file = stationary-outside +#@ args = method=realspace +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="1.2" +#@ args = max_step=4 +#@ +#@ [3d-spectral-linear] +#@ dim = 3 +#@ check-file = linear +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" +#@ args = max_step=1 +#@ +#@ [3d-spectral-linear-floor] +#@ dim = 3 +#@ check-file = linear-floor +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0125 +#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" +#@ args = max_step=1 +#@ +#@ [3d-spectral-nonlinear] +#@ dim = 3 +#@ check-file = nonlinear +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = max_step=40 +#@ +#@ [3d-spectral-stationary] +#@ dim = 3 +#@ check-file = stationary +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="-1" +#@ args = max_step=4 +#@ +#@ [3d-spectral-stationary-outside] +#@ dim = 3 +#@ check-file = stationary-outside +#@ fft = true +#@ args = method=spectral +#@ args = mobility_floor=0.0 +#@ args = eta.ic.expression.region0="1.2" +#@ args = max_step=4 +#@ +#@ [2d-realspace-mpi] +#@ dim = 2 +#@ nprocs = 2 +#@ check-file = nonlinear +#@ args = method=realspace +#@ args = max_step=40 +#@ +#@ [2d-spectral-mpi] +#@ dim = 2 +#@ nprocs = 2 +#@ check-file = nonlinear +#@ args = method=spectral +#@ args = max_step=40 +#@ fft = true +#@ +#@ [2d-spectral-amr] +#@ dim = 2 +#@ fft = true +#@ check-file = amr +#@ args = method=spectral amr.max_level=1 amr.nsubsteps=1 +#@ args = refinement_threshold=0.015 amr.regrid_int=2 max_step=8 +#@ args = eta.ic.expression.region0="0.15+0.1*exp(-40*((x-0.5)^2+(y-0.75)^2))" + +alamo.program = cahnhilliard +plot_file = tests/SinglyDegenerateCahnHilliard/output +amr.max_level = 0 +amr.n_cell = 32 32 32 +amr.max_grid_size = 16 +amr.blocking_factor = 4 +amr.plot_int = 1 +geometry.prob_lo = 0 0 0 +geometry.prob_hi = 1 1.5 2 +geometry.is_periodic = 1 1 1 +eta.bc.constant.type.xlo = periodic +eta.bc.constant.type.xhi = periodic +eta.bc.constant.type.ylo = periodic +eta.bc.constant.type.yhi = periodic +eta.bc.constant.type.zlo = periodic +eta.bc.constant.type.zhi = periodic +eta.ic.type = expression +eta.ic.expression.constant.pi = 3.14159265358979323846 +eta.ic.expression.region0 = "0.15+0.1*cos(2*pi*x)*cos(2*pi*y/1.5)" +method = realspace +mobility = singly_degenerate +mobility_floor = 0 +L = 1 +gamma = 0.01 +timestep = 0.00001 +max_step = 40 +stop_time = 1 diff --git a/tests/SinglyDegenerateCahnHilliard/test b/tests/SinglyDegenerateCahnHilliard/test new file mode 100755 index 0000000000..f63daf2963 --- /dev/null +++ b/tests/SinglyDegenerateCahnHilliard/test @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +"""Check conservation, dissipation, pure phases, and a linearized Fourier mode.""" +import sys +from pathlib import Path + +import numpy as np +import yt + +yt.set_log_level(50) +outdir = Path(sys.argv[1]) +case = sys.argv[2] +plots = sorted(outdir.glob("*cell")) +assert len(plots) >= 2, "Missing initial/final plotfiles" + + +def read(path): + ds = yt.load(str(path)) + level = ds.index.max_level + dimensions = ds.domain_dimensions.copy() + dimensions[:ds.dimensionality] *= 2**level + grid = ds.covering_grid(level, ds.domain_left_edge, dimensions) + eta = np.asarray(grid[("boxlib", "eta")]).squeeze() + dx = np.asarray(ds.domain_width)[:eta.ndim] / np.array(eta.shape) + return ds, eta, dx + + +def energy(eta, dx): + density = 0.25 * (eta**2 - 1)**2 + for axis, h in enumerate(dx): + density += 0.005 * ((np.roll(eta, -1, axis) - eta) / h)**2 + return density.mean() + + +def mean(ds): + cells = ds.all_data() + volumes = cells[("index", "cell_volume")] + return float((cells[("boxlib", "eta")] * volumes).sum() / volumes.sum()) + + +initial_ds, initial, dx = read(plots[0]) +final_ds, final, _ = read(plots[-1]) +assert np.isfinite(final).all(), "Nonfinite SDCH solution" +if case != "amr-real": + assert abs(mean(final_ds) - mean(initial_ds)) < 2e-12, "SDCH changed the mean" +assert float(final_ds.current_time) > float(initial_ds.current_time) + +if case.startswith("stationary"): + np.testing.assert_allclose(final, initial, rtol=0, atol=2e-13) +elif case.startswith("linear"): + # Linearize around eta=0.15. A product of cosines is a Laplacian eigenmode. + base = 0.15 + alpha = (base + 1) / 2 + floor = 0.0125 if case.endswith("floor") else 0.0 + mobility = alpha**2 * (1-alpha)**2 + floor + widths = np.asarray(initial_ds.domain_width)[:initial.ndim] + modes = np.ones(initial.ndim) + if case == "linear-odd": + modes[1] = 15 + wave = 2*np.pi*modes/widths + lap = np.sum(4*np.sin(wave*dx/2)**2/dx**2) + curvature = 3*base**2 - 1 + 0.01*lap + potential = np.asarray(final_ds.covering_grid(0, final_ds.domain_left_edge, + final_ds.domain_dimensions)[("boxlib", "int")]).squeeze() + expected_potential = base**3 - base + (initial-base) * curvature + np.testing.assert_allclose(potential, expected_potential, rtol=0, atol=1e-12) + spectral = "spectral" in outdir.name + if spectral: + stabilizer = 0.01 * (1/16 + floor) + rate = -mobility * np.sum(wave*np.sin(wave*dx)/dx) * curvature + rate /= 1 + 1e-5 * stabilizer * np.sum(wave**2)**2 + else: + rate = -mobility * lap * curvature + expected = base + (initial-base) * (1 + 1e-5*rate) + np.testing.assert_allclose(final, expected, rtol=0, atol=1e-13) + assert np.max(abs(final-initial)) > 1e-11, "Fourier mode did not evolve" +else: + assert np.max(abs(final-initial)) > 1e-7, "SDCH field did not evolve" + if case.startswith("amr"): + assert initial_ds.index.max_level == 1, "AMR regression did not refine" + fine_cells = sum(np.prod(g.ActiveDimensions) for g in initial_ds.index.grids if g.Level == 1) + assert 0 < fine_cells < np.prod(initial_ds.domain_dimensions[:initial.ndim] * 2), "AMR regression missed partial coverage" + else: + assert energy(final, dx) < energy(initial, dx), "SDCH failed to dissipate energy" +print(f"{case}: verified SDCH evolution") From a3f2f2fd58105e59a62d4a086475a65b7f5a75f6 Mon Sep 17 00:00:00 2001 From: Matthew Mehrtens <12023414+mcmehrtens@users.noreply.github.com> Date: Wed, 30 Sep 2026 13:20:50 -0500 Subject: [PATCH 2/4] Fix SDCH continuation indentation --- src/Integrator/CahnHilliard.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Integrator/CahnHilliard.cpp b/src/Integrator/CahnHilliard.cpp index abe31b6dae..2b3f8508f9 100644 --- a/src/Integrator/CahnHilliard.cpp +++ b/src/Integrator/CahnHilliard.cpp @@ -35,7 +35,7 @@ void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp) pp.query_default("mobility_floor", value.mobility_floor, 0.0); // Coefficient of the spectral fourth-order stabilizer. pp.query_default("spectral_stabilization", value.spectral_stabilization, - value.gamma * (value.L / 16.0 + value.mobility_floor)); + value.gamma * (value.L / 16.0 + value.mobility_floor)); // Regridding criterion pp.query_default("refinement_threshold",value.refinement_threshold, 1E100); @@ -131,7 +131,7 @@ Set::Scalar DegenerateMobility(Set::Scalar eta, Set::Scalar L, Set::Scalar floor } void ChemicalPotential(const amrex::MultiFab& eta_mf, amrex::MultiFab& mu_mf, - const amrex::Geometry& geometry, Set::Scalar gamma) + const amrex::Geometry& geometry, Set::Scalar gamma) { const auto dx = geometry.CellSizeArray(); for (amrex::MFIter mfi(eta_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) @@ -249,8 +249,8 @@ CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar if (length[2] % 2 == 0 && p == length[2]/2) k3 = 0.0; #endif Set::Complex div_flux = AMREX_D_TERM(I * k1 * flux_hat(m,n,p,0), - + I * k2 * flux_hat(m,n,p,1), - + I * k3 * flux_hat(m,n,p,2)); + + I * k2 * flux_hat(m,n,p,1), + + I * k3 * flux_hat(m,n,p,2)); // Keeping the zero mode unchanged conserves the mean without clipping. eta_hat(m,n,p) += dt * div_flux / (1.0 + dt * S * omega2 * omega2); }); From a65bb0b82801f012ca4ba6707f67ea0adebebd61 Mon Sep 17 00:00:00 2001 From: Matthew Mehrtens <12023414+mcmehrtens@users.noreply.github.com> Date: Mon, 5 Oct 2026 13:17:10 -0500 Subject: [PATCH 3/4] Use compatible conservative SDCH operators and test FFT dimensions --- .github/workflows/linux.yml | 22 +++--- src/Integrator/CahnHilliard.H | 6 +- src/Integrator/CahnHilliard.cpp | 75 ++++++++++--------- src/Operator/Spectral/FFT.H | 26 ------- tests/SinglyDegenerateCahnHilliard/Readme.rst | 15 +++- tests/SinglyDegenerateCahnHilliard/input | 71 ++++++++++++++++++ tests/SinglyDegenerateCahnHilliard/test | 41 ++++++++-- 7 files changed, 175 insertions(+), 81 deletions(-) diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 69449f6ed9..b6b842d3ac 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -234,8 +234,12 @@ jobs: echo "## 📄 Report: test-2d-24-04-clang" >> $GITHUB_STEP_SUMMARY echo "" >> $GITHUB_STEP_SUMMARY echo "[⬇️ Click here to download](${{ steps.upload_report.outputs.artifact-url }})" >> $GITHUB_STEP_SUMMARY - test-2d-24-04-clang-fft: + test-24-04-clang-fft: + name: test-${{ matrix.dim }}d-24-04-clang-fft runs-on: ubuntu-24.04 + strategy: + matrix: + dim: [2, 3] steps: - uses: actions/checkout@v4 - name: installing preliminaries @@ -244,19 +248,19 @@ jobs: run: sudo apt install clang - name: installing fftw3 run: sudo apt install libfftw3-dev - - name: configure (2d) - run: ./configure --dim=2 --comp=clang++ --fft + - name: configure (${{ matrix.dim }}d) + run: ./configure --dim=${{ matrix.dim }} --comp=clang++ --fft - name: make run: make - name: regression tests - run: scripts/runtests.py --serial --no-coverage --comp=clang++ --timeout=200 --no-backspace --fft-only + run: scripts/runtests.py --serial --no-coverage --dim=${{ matrix.dim }} --comp=clang++ --timeout=200 --no-backspace --fft-only - name: create report if: always() - run: python3 scripts/package_report.py --name test-2d-24-04-clang-fft + run: python3 scripts/package_report.py --name test-${{ matrix.dim }}d-24-04-clang-fft - name: report artifact summary if: always() run: | - name="test-2d-24-04-clang-fft" + name="test-${{ matrix.dim }}d-24-04-clang-fft" echo "### Report bundle: $name" >> $GITHUB_STEP_SUMMARY if [ -d "$name" ]; then echo "Files: $(find "$name" -type f | wc -l)" >> $GITHUB_STEP_SUMMARY @@ -270,11 +274,11 @@ jobs: id: upload_report uses: actions/upload-artifact@v4 with: - name: test-2d-24-04-clang-fft - path: test-2d-24-04-clang-fft + name: test-${{ matrix.dim }}d-24-04-clang-fft + path: test-${{ matrix.dim }}d-24-04-clang-fft - name: Add report link to summary if: always() run: | - echo "## 📄 Report: test-2d-24-04-clang-fft" >> $GITHUB_STEP_SUMMARY + echo "## 📄 Report: test-${{ matrix.dim }}d-24-04-clang-fft" >> $GITHUB_STEP_SUMMARY echo "" >> $GITHUB_STEP_SUMMARY echo "[⬇️ Click here to download](${{ steps.upload_report.outputs.artifact-url }})" >> $GITHUB_STEP_SUMMARY diff --git a/src/Integrator/CahnHilliard.H b/src/Integrator/CahnHilliard.H index ad2a6897d0..37ab65890a 100644 --- a/src/Integrator/CahnHilliard.H +++ b/src/Integrator/CahnHilliard.H @@ -26,8 +26,10 @@ // // where :math:`\mu=\eta^3-\eta-\gamma\nabla^2\eta` and // :math:`M_0` is :code:`mobility_floor` (default zero). -// The realspace method uses conservative face fluxes; the spectral method -// stabilizes the flux divergence with :code:`spectral_stabilization`, default +// Both SDCH methods use conservative face fluxes and the same discrete Laplacian. +// The spectral method applies the Fourier multiplier :math:`1/(1+\Delta t S D^2)` +// to the increment, where :math:`D` is the negative discrete Laplacian symbol +// and :math:`S` is :code:`spectral_stabilization`, default // :math:`\gamma(L/16+M_0)`. SDCH requires periodic boundaries. Its spectral // method assembles a uniform finest-level grid and restricts the update to // the AMR hierarchy; realspace AMR has no flux refluxing at coarse/fine interfaces. diff --git a/src/Integrator/CahnHilliard.cpp b/src/Integrator/CahnHilliard.cpp index 2b3f8508f9..905eeda6ee 100644 --- a/src/Integrator/CahnHilliard.cpp +++ b/src/Integrator/CahnHilliard.cpp @@ -130,6 +130,25 @@ Set::Scalar DegenerateMobility(Set::Scalar eta, Set::Scalar L, Set::Scalar floor return L * alpha * alpha * (1.0 - alpha) * (1.0 - alpha) + floor; } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +Set::Scalar ConservativeRHS(const amrex::Array4& eta, + const amrex::Array4& mu, + int i, int j, int k, const Set::Scalar* dx, + Set::Scalar L, Set::Scalar floor) +{ + const Set::Scalar M = DegenerateMobility(eta(i,j,k), L, floor); + Set::Scalar rhs = 0.0; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) + { + const int di = dir == 0, dj = dir == 1, dk = dir == 2; + const Set::Scalar plus = 0.5 * (M + DegenerateMobility(eta(i+di,j+dj,k+dk), L, floor)); + const Set::Scalar minus = 0.5 * (M + DegenerateMobility(eta(i-di,j-dj,k-dk), L, floor)); + rhs += (plus * (mu(i+di,j+dj,k+dk) - mu(i,j,k)) + - minus * (mu(i,j,k) - mu(i-di,j-dj,k-dk))) / (dx[dir] * dx[dir]); + } + return rhs; +} + void ChemicalPotential(const amrex::MultiFab& eta_mf, amrex::MultiFab& mu_mf, const amrex::Geometry& geometry, Set::Scalar gamma) { @@ -176,17 +195,7 @@ CahnHilliard::AdvanceDegenerateReal (int lev, Set::Scalar time, Set::Scalar dt) const auto next = etanew_mf[lev]->array(mfi); amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) { - const Set::Scalar M = DegenerateMobility(eta(i,j,k), L_local, floor); - Set::Scalar rhs = 0.0; - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - { - const int di = dir == 0, dj = dir == 1, dk = dir == 2; - const Set::Scalar plus = 0.5 * (M + DegenerateMobility(eta(i+di,j+dj,k+dk), L_local, floor)); - const Set::Scalar minus = 0.5 * (M + DegenerateMobility(eta(i-di,j-dj,k-dk), L_local, floor)); - rhs += (plus * (mu(i+di,j+dj,k+dk) - mu(i,j,k)) - - minus * (mu(i,j,k) - mu(i-di,j-dj,k-dk))) / (dx[dir] * dx[dir]); - } - next(i,j,k) = eta(i,j,k) + dt * rhs; + next(i,j,k) = eta(i,j,k) + dt * ConservativeRHS(eta, mu, i, j, k, dx.data(), L_local, floor); }); } etanew_mf[lev]->FillBoundary(geom[lev].periodicity()); @@ -209,7 +218,7 @@ CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar amrex::Vector potential_hierarchy(lev + 1); for (int ilev = 0; ilev <= lev; ++ilev) potential_hierarchy[ilev] = intermediate[ilev].get(); FFT::UniformToComposite(mu_mf, potential_hierarchy, geom, refRatio(), lev, 1); - amrex::MultiFab flux_mf(eta_mf->boxArray(), eta_mf->DistributionMap(), AMREX_SPACEDIM, 0); + amrex::MultiFab rhs_mf(eta_mf->boxArray(), eta_mf->DistributionMap(), 1, 0); const auto dx = geom[lev].CellSizeArray(); const Set::Scalar L_local = L; const Set::Scalar floor = mobility_floor; @@ -217,42 +226,38 @@ CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar { const auto eta = eta_mf->const_array(mfi); const auto mu = mu_mf.const_array(mfi); - const auto flux = flux_mf.array(mfi); + const auto rhs = rhs_mf.array(mfi); amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) { - Set::Scalar M = DegenerateMobility(eta(i,j,k), L_local, floor); - Set::Vector grad_mu = Numeric::Gradient(mu, i, j, k, 0, dx.data()); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) flux(i,j,k,dir) = M * grad_mu(dir); + // Face gradients and their adjoint divergence dissipate the discrete energy. + rhs(i,j,k) = ConservativeRHS(eta, mu, i, j, k, dx.data(), L_local, floor); }); } auto eta_hat_mf = fft.MakeSpectralFab(); - auto flux_hat_mf = fft.MakeSpectralFab(AMREX_SPACEDIM); + auto rhs_hat_mf = fft.MakeSpectralFab(); fft.Forward(*eta_mf, eta_hat_mf); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) fft.Forward(flux_mf, flux_hat_mf, dir, dir); + fft.Forward(rhs_mf, rhs_hat_mf); const Set::Scalar S = spectral_stabilization; - const Set::Complex I(0.0, 1.0); const auto length = geom[lev].Domain().length3d(); for (amrex::MFIter mfi(eta_hat_mf, false); mfi.isValid(); ++mfi) { const auto eta_hat = eta_hat_mf.array(mfi); - const auto flux_hat = flux_hat_mf.const_array(mfi); - fft.ParallelForModes(mfi.tilebox(), [=] AMREX_GPU_DEVICE( - int m, int n, int p, Set::Scalar k1, Set::Scalar k2, Set::Scalar k3, Set::Scalar omega2) + const auto rhs_hat = rhs_hat_mf.const_array(mfi); + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int m, int n, int p) { - amrex::ignore_unused(k2, k3); - // A real Fourier derivative vanishes at even-grid Nyquist modes. - if (length[0] % 2 == 0 && m == length[0]/2) k1 = 0.0; -#if AMREX_SPACEDIM >= 2 - if (length[1] % 2 == 0 && n == length[1]/2) k2 = 0.0; -#endif -#if AMREX_SPACEDIM >= 3 - if (length[2] % 2 == 0 && p == length[2]/2) k3 = 0.0; + // The face-flux Laplacian symbol includes both odd-grid and Nyquist modes. + const Set::Scalar sx = std::sin(Set::Constant::Pi * m / length[0]); + const Set::Scalar sy = std::sin(Set::Constant::Pi * n / length[1]); + amrex::ignore_unused(p); +#if AMREX_SPACEDIM == 3 + const Set::Scalar sz = std::sin(Set::Constant::Pi * p / length[2]); #endif - Set::Complex div_flux = AMREX_D_TERM(I * k1 * flux_hat(m,n,p,0), - + I * k2 * flux_hat(m,n,p,1), - + I * k3 * flux_hat(m,n,p,2)); - // Keeping the zero mode unchanged conserves the mean without clipping. - eta_hat(m,n,p) += dt * div_flux / (1.0 + dt * S * omega2 * omega2); + const Set::Scalar lap = AMREX_D_TERM(4.0 * sx*sx / (dx[0]*dx[0]), + + 4.0 * sy*sy / (dx[1]*dx[1]), + + 4.0 * sz*sz / (dx[2]*dx[2])); + // Set the conserved mode exactly; the face flux sum is zero up to roundoff. + if (m != 0 || n != 0 || p != 0) + eta_hat(m,n,p) += dt * rhs_hat(m,n,p) / (1.0 + dt * S * lap * lap); }); } fft.Backward(eta_hat_mf, *eta_mf); diff --git a/src/Operator/Spectral/FFT.H b/src/Operator/Spectral/FFT.H index dafa0546d9..c7210f70af 100644 --- a/src/Operator/Spectral/FFT.H +++ b/src/Operator/Spectral/FFT.H @@ -119,32 +119,6 @@ public: }); } - template - void ParallelForModes (const amrex::Box& bx, F&& f) const - { - int length[AMREX_SPACEDIM]; - Set::Scalar kscale[AMREX_SPACEDIM]; - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - { - length[dir] = m_length[dir]; - kscale[dir] = m_kscale[dir]; - } - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int m, int n, int p) { - Set::Scalar k1 = K(0, m, length, kscale); - Set::Scalar k2 = 0.0; - Set::Scalar k3 = 0.0; -#if AMREX_SPACEDIM >= 2 - k2 = (n <= (length[1]-1)/2 ? n : n-length[1]) * kscale[1]; -#endif -#if AMREX_SPACEDIM >= 3 - k3 = (p <= (length[2]-1)/2 ? p : p-length[2]) * kscale[2]; -#endif - Set::Scalar omega2 = k1*k1 + k2*k2 + k3*k3; - f(m, n, p, k1, k2, k3, omega2); - }); - } - void Forward (const amrex::MultiFab& in, amrex::FabArray >& out, int scomp = 0, int dcomp = 0) { m_fft.forward(in, out, scomp, dcomp); diff --git a/tests/SinglyDegenerateCahnHilliard/Readme.rst b/tests/SinglyDegenerateCahnHilliard/Readme.rst index 645f1dbbbe..766d9d5f8a 100644 --- a/tests/SinglyDegenerateCahnHilliard/Readme.rst +++ b/tests/SinglyDegenerateCahnHilliard/Readme.rst @@ -4,13 +4,20 @@ Use :code:`mobility = singly_degenerate` with :code:`method = realspace` or :code:`method = spectral` (requires a build configured with :code:`--fft`). The default :code:`mobility_floor = 0` makes the pure phases immobile; :code:`spectral_stabilization` defaults to :math:`\gamma(L/16+M_0)`. -Constant mobility remains the default. Periodic boundary conditions are required. The realspace timestep must satisfy -the explicit fourth-order stability constraint; neither method clips the solution. +Both methods use face fluxes with arithmetic face mobility and a centered +Laplacian for the chemical potential. The spectral method uses FFTs to apply +a discrete biharmonic stabilizer to the increment, including Nyquist modes. +It does not guarantee energy decrease for arbitrary timesteps. +Constant mobility remains the default. Periodic boundary conditions are required. +The realspace timestep must satisfy the explicit fourth-order stability constraint; +neither method clips the solution. The regression sections check a linearized Fourier mode on anisotropic domains, nonzero mobility floors, stationary phases, conservation, and energy dissipation -in two and three dimensions. MPI sections exercise patch/rank boundaries, and -the spectral AMR section exercises a partially refined grid with regridding. +in two and three dimensions. Grid-scale checkerboards must decay, and a +strongly varying mobility case must dissipate energy at small timesteps. MPI +sections exercise patch/rank boundaries, and the spectral AMR section exercises +a partially refined grid with regridding. Realspace AMR does not reflux coarse/fine fluxes; use a uniform grid when exact mass conservation is required. Spectral AMR assembles a uniform finest-level grid, which requires memory proportional to that entire grid. diff --git a/tests/SinglyDegenerateCahnHilliard/input b/tests/SinglyDegenerateCahnHilliard/input index 6f3a3a95b5..f32de6f8bf 100644 --- a/tests/SinglyDegenerateCahnHilliard/input +++ b/tests/SinglyDegenerateCahnHilliard/input @@ -207,6 +207,77 @@ #@ args = refinement_threshold=0.015 amr.regrid_int=2 max_step=8 #@ args = eta.ic.expression.region0="0.15+0.1*exp(-40*((x-0.5)^2+(y-0.75)^2))" +#@ [2d-spectral-checkerboard-x] +#@ dim = 2 +#@ fft = true +#@ check-file = checkerboard-x +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)" +#@ +#@ [2d-spectral-checkerboard-y] +#@ dim = 2 +#@ fft = true +#@ check-file = checkerboard-y +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*y/1.5)" +#@ +#@ [2d-spectral-checkerboard-xy] +#@ dim = 2 +#@ fft = true +#@ check-file = checkerboard-xy +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)*sin(32*pi*y/1.5)" +#@ +#@ [2d-spectral-strong-mobility] +#@ dim = 2 +#@ args = timestep=1e-8 +#@ +#@ [2d-spectral-strong-mobility-small] +#@ dim = 2 +#@ args = timestep=1e-10 +#@ +#@ [3d-spectral-checkerboard-x] +#@ dim = 3 +#@ fft = true +#@ check-file = checkerboard-x +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)" +#@ +#@ [3d-spectral-checkerboard-y] +#@ dim = 3 +#@ fft = true +#@ check-file = checkerboard-y +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*y/1.5)" +#@ +#@ [3d-spectral-checkerboard-z] +#@ dim = 3 +#@ fft = true +#@ check-file = checkerboard-z +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*z/2)" +#@ +#@ [3d-spectral-checkerboard-xyz] +#@ dim = 3 +#@ fft = true +#@ check-file = checkerboard-xyz +#@ args = method=spectral max_step=4 +#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)*sin(32*pi*y/1.5)*sin(32*pi*z/2)" +#@ +#@ [3d-spectral-strong-mobility] +#@ dim = 3 +#@ args = timestep=1e-8 +#@ +#@ [3d-spectral-strong-mobility-small] +#@ dim = 3 +#@ args = timestep=1e-10 +#@ +#@ [*spectral-strong-mobility*] +#@ fft = true +#@ check-file = strong-mobility +#@ args = method=spectral amr.n_cell=16 16 16 max_step=1 +#@ args = eta.ic.expression.region0="(0.17623171536106855)*(x>=0 and x<0.0625)+(0.9999915189178985)*(x>=0.0625 and x<0.125)+(0.99957222732831352)*(x>=0.125 and x<0.1875)+(-0.33828120309429632)*(x>=0.1875 and x<0.25)+(0.96308521781430201)*(x>=0.25 and x<0.3125)+(0.86409632475355058)*(x>=0.3125 and x<0.375)+(0.98990644052604582)*(x>=0.375 and x<0.4375)+(-0.98167620233468933)*(x>=0.4375 and x<0.5)+(0.99984268465982185)*(x>=0.5 and x<0.5625)+(0.99898999623396001)*(x>=0.5625 and x<0.625)+(-0.54394735059310761)*(x>=0.625 and x<0.6875)+(-0.99762312075226711)*(x>=0.6875 and x<0.75)+(-0.99954221397642695)*(x>=0.75 and x<0.8125)+(0.99930823691403747)*(x>=0.8125 and x<0.875)+(0.99991743650726583)*(x>=0.875 and x<0.9375)+(0.89293086073486305)*(x>=0.9375 and x<1)" + alamo.program = cahnhilliard plot_file = tests/SinglyDegenerateCahnHilliard/output amr.max_level = 0 diff --git a/tests/SinglyDegenerateCahnHilliard/test b/tests/SinglyDegenerateCahnHilliard/test index f63daf2963..b12770a396 100755 --- a/tests/SinglyDegenerateCahnHilliard/test +++ b/tests/SinglyDegenerateCahnHilliard/test @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -"""Check conservation, dissipation, pure phases, and a linearized Fourier mode.""" +"""Check conservation, dissipation, pure phases, and Fourier-mode evolution.""" import sys from pathlib import Path @@ -64,15 +64,46 @@ elif case.startswith("linear"): expected_potential = base**3 - base + (initial-base) * curvature np.testing.assert_allclose(potential, expected_potential, rtol=0, atol=1e-12) spectral = "spectral" in outdir.name + rate = -mobility * lap * curvature if spectral: stabilizer = 0.01 * (1/16 + floor) - rate = -mobility * np.sum(wave*np.sin(wave*dx)/dx) * curvature - rate /= 1 + 1e-5 * stabilizer * np.sum(wave**2)**2 - else: - rate = -mobility * lap * curvature + rate /= 1 + 1e-5 * stabilizer * lap**2 expected = base + (initial-base) * (1 + 1e-5*rate) np.testing.assert_allclose(final, expected, rtol=0, atol=1e-13) assert np.max(abs(final-initial)) > 1e-11, "Fourier mode did not evolve" +elif case == "strong-mobility": + # This grid-scale field previously evolved uphill even as dt tended to zero. + mobility = ((1-initial**2) / 4)**2 + assert mobility.max() / mobility.min() > 1e6, "Weak mobility variation" + mu = initial**3 - initial + for axis, h in enumerate(dx): + mu -= 0.01 * (np.roll(initial, -1, axis) - 2*initial + + np.roll(initial, 1, axis)) / h**2 + dissipation = 0.0 + for axis, h in enumerate(dx): + face_mobility = (mobility + np.roll(mobility, -1, axis)) / 2 + dissipation += np.mean(face_mobility * ((np.roll(mu, -1, axis)-mu)/h)**2) + dt = float(final_ds.current_time - initial_ds.current_time) + energies = [energy(read(path)[1], dx) for path in plots] + assert np.all(np.diff(energies) < 0), "Strong-mobility SDCH increased energy" + observed = (energies[-1] - energies[0]) / dt + # At small dt, the discrete energy derivative is minus the squared face gradients. + np.testing.assert_allclose(observed, -dissipation, rtol=2e-3, atol=0) + print(f"energy change: {energies[-1]-energies[0]:.12g}; rate: {observed:.12g}") +elif case.startswith("checkerboard"): + # For eta=b+a*s, s=+/-1, every active face sees the same mean mobility. + base, amplitude = 0.15, 0.01 + axes = case.split("-")[-1] + lap = sum(4 / dx["xyz".index(axis)]**2 for axis in axes) + dt = 1e-5 + steps = round(float(final_ds.current_time) / dt) + for _ in range(steps): + mobility = sum(((1-eta**2)/4)**2 for eta in (base+amplitude, base-amplitude)) / 2 + curvature = 3*base**2 - 1 + amplitude**2 + 0.01*lap + amplitude *= 1 - dt*mobility*lap*curvature / (1 + dt*0.01/16*lap**2) + np.testing.assert_allclose(final, base+(initial-base)*amplitude/0.01, rtol=0, atol=2e-13) + assert np.max(abs(final-base)) < 0.0099, "Checkerboard mode did not decay" + assert energy(final, dx) < energy(initial, dx), "Checkerboard energy did not decay" else: assert np.max(abs(final-initial)) > 1e-7, "SDCH field did not evolve" if case.startswith("amr"): From 7c6416304b04276b0ca37a6d6bcde31cf6e55e95 Mon Sep 17 00:00:00 2001 From: Matthew Mehrtens <12023414+mcmehrtens@users.noreply.github.com> Date: Mon, 26 Oct 2026 15:10:11 -0500 Subject: [PATCH 4/4] Port reusable Cahn-Hilliard updates and focus regressions --- scripts/runtests.py | 2 + src/Integrator/CahnHilliard.H | 37 ++- src/Integrator/CahnHilliard.cpp | 233 ++++++++++++++- tests/SinglyDegenerateCahnHilliard/Readme.rst | 26 +- tests/SinglyDegenerateCahnHilliard/input | 273 +++--------------- tests/SinglyDegenerateCahnHilliard/test | 142 ++++++++- 6 files changed, 443 insertions(+), 270 deletions(-) diff --git a/scripts/runtests.py b/scripts/runtests.py index 23449f8cec..80488a3d25 100755 --- a/scripts/runtests.py +++ b/scripts/runtests.py @@ -145,6 +145,7 @@ def readMetadata(metadatafile): parser.add_argument('--mpirun-flags',dest="mpirun_flags",default="",help="Extra arguments to pass to mpirun (like --oversubscribe). All arguments must be in a string.") parser.add_argument('--fft',dest="fft",default=False,action='store_true',help="Enable fft-based tests") parser.add_argument('--fft-only',dest="fft_only",default=False,action='store_true',help="Run fft tests only") +parser.add_argument('--hdf5', default=False, action='store_true', help='Use HDF5-enabled executables') parser.add_argument('--cuda', const='auto', default=None, nargs='?', type=str, help="Run CUDA tests, optionally specifying the executable architecture suffix") parser.add_argument('--post-timeout', dest="post_timeout", default=10000, help='How long to wait before skipping results posting') @@ -445,6 +446,7 @@ def test(testdir): if args.profile: exestr += "-profile" if args.perf: exestr += "-perf" if coverage: exestr += "-coverage" + if args.hdf5: exestr += "-hdf5" if args.cuda: cuda_arch = args.cuda if cuda_arch == "auto": diff --git a/src/Integrator/CahnHilliard.H b/src/Integrator/CahnHilliard.H index 37ab65890a..aea144cc49 100644 --- a/src/Integrator/CahnHilliard.H +++ b/src/Integrator/CahnHilliard.H @@ -33,7 +33,15 @@ // :math:`\gamma(L/16+M_0)`. SDCH requires periodic boundaries. Its spectral // method assembles a uniform finest-level grid and restricts the update to // the AMR hierarchy; realspace AMR has no flux refluxing at coarse/fine interfaces. -// The solution is not clipped, so timesteps must resolve the interface dynamics. +// By default the solution is not clipped, so timesteps must resolve the +// interface dynamics. Clipping can change the conserved mean. +// Optional input and field names are :code:`input_name`, :code:`field_name`, +// :code:`old_field_name`, and :code:`intermediate_name`. SDCH can optionally +// clip its updated field with :code:`clipping` (default false). Derived integrators can override +// :code:`FillMobilityMask` to apply a spatial multiplier to the face mobility. +// The hook receives each level's native box layout and valid cells; the caller +// fills periodic ghosts. Set :code:`diagnostics = true` for SDCH mobility, mask, flux, RHS, raw update, +// and update delta fields alongside the solution. // // This is tested in :ref:`CahnHilliard` and :ref:`SinglyDegenerateCahnHilliard` @@ -51,7 +59,7 @@ namespace Integrator { -class CahnHilliard : public Integrator +class CahnHilliard : virtual public Integrator { public: static constexpr const char* name = "cahnhilliard"; @@ -86,14 +94,25 @@ protected: /// Mark any cells that need to be refined void TagCellsForRefinement (int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow) override; -private: + /// Fill a cell-centered multiplier for the singly-degenerate mobility. + /// The default is one everywhere; derived integrators may override it. + virtual void FillMobilityMask(int lev, amrex::MultiFab& mask); + void InitializeDiagnostics(int lev); + +protected: Set::Field etanew_mf; /// The new value for eta this timestep Set::Field etaold_mf; /// Last timestep's value for eta Set::Field intermediate; /// Intermediate field used for CH kinetics + Set::Field diagnostic_mobility; + Set::Field diagnostic_mobility_mask; + Set::Field diagnostic_flux; + Set::Field diagnostic_rhs; + Set::Field diagnostic_raw; + Set::Field diagnostic_delta; - BC::BC *bc; /// eta's bc object - IC::IC *ic; /// eta's ic object + BC::BC *bc = nullptr; /// eta's bc object + IC::IC *ic = nullptr; /// eta's ic object Set::Scalar gamma = NAN; Set::Scalar L = NAN; @@ -103,6 +122,14 @@ private: std::string mobility = "constant"; Set::Scalar mobility_floor = 0.0; Set::Scalar spectral_stabilization = NAN; + Set::Scalar tstart = 1E100; + Set::Scalar tstart_timestep = -1.0; + bool clipping = false; + bool output_diagnostics = false; + std::string input_name = "eta"; + std::string field_name = "eta"; + std::string old_field_name = "eta_old"; + std::string intermediate_name = "int"; }; } #endif diff --git a/src/Integrator/CahnHilliard.cpp b/src/Integrator/CahnHilliard.cpp index 905eeda6ee..4d70743847 100644 --- a/src/Integrator/CahnHilliard.cpp +++ b/src/Integrator/CahnHilliard.cpp @@ -1,6 +1,8 @@ #include #include #include +#include +#include #include "IC/Expression.H" @@ -39,10 +41,19 @@ void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp) // Regridding criterion pp.query_default("refinement_threshold",value.refinement_threshold, 1E100); + pp.query_default("input_name", value.input_name, value.input_name); + pp.query_default("field_name", value.field_name, value.field_name); + pp.query_default("old_field_name", value.old_field_name, value.old_field_name); + pp.query_default("intermediate_name", value.intermediate_name, value.intermediate_name); + pp.query_default("clipping", value.clipping, false); + pp.query_default("diagnostics", value.output_diagnostics, false); + pp.query_default("tstart", value.tstart, 1E100); + pp.query_default("tstart_timestep", value.tstart_timestep, -1.0); + // initial condition for :math:`\eta` - pp.select_default("eta.ic", value.ic, pp.forward_args(value.geom)); + pp.select_default(value.input_name + ".ic", value.ic, pp.forward_args(value.geom)); // boundary condition for :math:`\eta` - pp.select_default("eta.bc", value.bc, pp.forward_args(1)); + pp.select_default(value.input_name + ".bc", value.bc, pp.forward_args(1)); // Which method to use - realspace or spectral method. pp.query_validate("method",value.method,{"realspace","spectral"}); @@ -57,17 +68,47 @@ void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp) if (!value.geom[0].isAllPeriodic()) Util::ParmParseException(INFO, "mobility", "SDCH requires periodic boundaries"); } + if (value.output_diagnostics && value.mobility != "singly_degenerate" && !pp.InTraversalMode()) + Util::ParmParseException(INFO, "diagnostics", "CahnHilliard diagnostics require mobility=singly_degenerate"); - value.RegisterNewFab(value.etanew_mf, value.bc, 1, 1, "eta",true); - value.RegisterNewFab(value.intermediate, value.bc, 1, 1, "int",true); + value.RegisterNewFab(value.etanew_mf, value.bc, 1, 1, value.field_name,true); + value.RegisterNewFab(value.intermediate, value.bc, 1, 1, value.intermediate_name,true); + + if (value.output_diagnostics) + { + std::string prefix; + const auto separator = value.field_name.find_last_of('.'); + if (separator != std::string::npos) prefix = value.field_name.substr(0, separator + 1); + std::vector flux_suffix{"_x"}; +#if AMREX_SPACEDIM >= 2 + flux_suffix.push_back("_y"); +#endif +#if AMREX_SPACEDIM >= 3 + flux_suffix.push_back("_z"); +#endif + value.RegisterNewFab(value.diagnostic_mobility, nullptr, 1, 0, prefix + "mobility", true, false); + value.RegisterNewFab(value.diagnostic_mobility_mask, nullptr, 1, 0, prefix + "mobility_mask", true, false); + value.RegisterNewFab(value.diagnostic_flux, nullptr, AMREX_SPACEDIM, 0, prefix + "flux", true, false, flux_suffix); + value.RegisterNewFab(value.diagnostic_rhs, nullptr, 1, 0, prefix + "rhs", true, false); + value.RegisterNewFab(value.diagnostic_raw, nullptr, 1, 0, prefix + "alpha_raw", true, false); + value.RegisterNewFab(value.diagnostic_delta, nullptr, 1, 0, prefix + "alpha_delta", true, false); + } if (value.method == "realspace") - value.RegisterNewFab(value.etaold_mf, value.bc, 1, 1, "eta_old",false); + value.RegisterNewFab(value.etaold_mf, value.bc, 1, 1, value.old_field_name,false); +} + +void +CahnHilliard::FillMobilityMask(int /*lev*/, amrex::MultiFab& mask) +{ + mask.setVal(1.0); } void CahnHilliard::Advance(int lev, Set::Scalar time, Set::Scalar dt) { + if (method == "realspace" && time > tstart && tstart_timestep > 0.0) + SetTimestep(tstart_timestep); if (mobility == "singly_degenerate") { if (method == "realspace") @@ -130,19 +171,40 @@ Set::Scalar DegenerateMobility(Set::Scalar eta, Set::Scalar L, Set::Scalar floor return L * alpha * alpha * (1.0 - alpha) * (1.0 - alpha) + floor; } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void CellFlux(const amrex::Array4& eta, + const amrex::Array4& mu, + const amrex::Array4& mask, + int i, int j, int k, const Set::Scalar* dx, + Set::Scalar L, Set::Scalar floor, + const amrex::Array4& flux) +{ + const Set::Scalar M = mask(i,j,k) * DegenerateMobility(eta(i,j,k), L, floor); + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) + { + const int di = dir == 0, dj = dir == 1, dk = dir == 2; + const Set::Scalar Mp = mask(i+di,j+dj,k+dk) * DegenerateMobility(eta(i+di,j+dj,k+dk),L,floor); + const Set::Scalar Mm = mask(i-di,j-dj,k-dk) * DegenerateMobility(eta(i-di,j-dj,k-dk),L,floor); + const Set::Scalar fp = 0.5 * (M + Mp) * (mu(i+di,j+dj,k+dk)-mu(i,j,k)) / dx[dir]; + const Set::Scalar fm = 0.5 * (M + Mm) * (mu(i,j,k)-mu(i-di,j-dj,k-dk)) / dx[dir]; + flux(i,j,k,dir) = 0.5 * (fp + fm); + } +} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Set::Scalar ConservativeRHS(const amrex::Array4& eta, const amrex::Array4& mu, + const amrex::Array4& mask, int i, int j, int k, const Set::Scalar* dx, Set::Scalar L, Set::Scalar floor) { - const Set::Scalar M = DegenerateMobility(eta(i,j,k), L, floor); + const Set::Scalar M = mask(i,j,k) * DegenerateMobility(eta(i,j,k), L, floor); Set::Scalar rhs = 0.0; for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { const int di = dir == 0, dj = dir == 1, dk = dir == 2; - const Set::Scalar plus = 0.5 * (M + DegenerateMobility(eta(i+di,j+dj,k+dk), L, floor)); - const Set::Scalar minus = 0.5 * (M + DegenerateMobility(eta(i-di,j-dj,k-dk), L, floor)); + const Set::Scalar plus = 0.5 * (M + mask(i+di,j+dj,k+dk) * DegenerateMobility(eta(i+di,j+dj,k+dk), L, floor)); + const Set::Scalar minus = 0.5 * (M + mask(i-di,j-dj,k-dk) * DegenerateMobility(eta(i-di,j-dj,k-dk), L, floor)); rhs += (plus * (mu(i+di,j+dj,k+dk) - mu(i,j,k)) - minus * (mu(i,j,k) - mu(i-di,j-dj,k-dk))) / (dx[dir] * dx[dir]); } @@ -173,6 +235,19 @@ CahnHilliard::AdvanceDegenerateReal (int lev, Set::Scalar time, Set::Scalar dt) etaold_mf[lev]->FillBoundary(geom[lev].periodicity()); ChemicalPotential(*etaold_mf[lev], *intermediate[lev], geom[lev], gamma); intermediate[lev]->FillBoundary(geom[lev].periodicity()); + amrex::MultiFab mobility_mask(etanew_mf[lev]->boxArray(), etanew_mf[lev]->DistributionMap(), 1, 1); + mobility_mask.setVal(1.0); + FillMobilityMask(lev, mobility_mask); + mobility_mask.FillBoundary(geom[lev].periodicity()); + std::unique_ptr coarse_mobility_mask; + if (lev > 0) + { + coarse_mobility_mask = std::make_unique( + etanew_mf[lev-1]->boxArray(), etanew_mf[lev-1]->DistributionMap(), 1, 1); + coarse_mobility_mask->setVal(1.0); + FillMobilityMask(lev-1, *coarse_mobility_mask); + coarse_mobility_mask->FillBoundary(geom[lev-1].periodicity()); + } if (lev > 0) { amrex::Vector coarse{intermediate[lev-1].get()}; @@ -183,19 +258,48 @@ CahnHilliard::AdvanceDegenerateReal (int lev, Set::Scalar time, Set::Scalar dt) amrex::FillPatchTwoLevels(*intermediate[lev], time, coarse, times, fine, times, 0, 0, 1, geom[lev-1], geom[lev], *bc, 0, *bc, 0, refRatio(lev-1), &amrex::cell_cons_interp, bcs, 0); + + amrex::Vector coarse_mask{coarse_mobility_mask.get()}; + amrex::Vector fine_mask{&mobility_mask}; + amrex::FillPatchTwoLevels(mobility_mask, time, coarse_mask, times, fine_mask, times, + 0, 0, 1, geom[lev-1], geom[lev], *bc, 0, *bc, 0, refRatio(lev-1), + &amrex::cell_cons_interp, bcs, 0); } const auto dx = geom[lev].CellSizeArray(); const Set::Scalar L_local = L; const Set::Scalar floor = mobility_floor; + const bool clipping_local = clipping; + const bool diagnostics_local = output_diagnostics; for (amrex::MFIter mfi(*etanew_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { const auto eta = etaold_mf[lev]->const_array(mfi); const auto mu = intermediate[lev]->const_array(mfi); + const auto mask = mobility_mask.const_array(mfi); const auto next = etanew_mf[lev]->array(mfi); + amrex::Array4 out_m, out_mask, out_raw, out_delta, out_rhs, out_flux; + if (diagnostics_local) + { + out_m = diagnostic_mobility[lev]->array(mfi); + out_mask = diagnostic_mobility_mask[lev]->array(mfi); + out_raw = diagnostic_raw[lev]->array(mfi); + out_delta = diagnostic_delta[lev]->array(mfi); + out_rhs = diagnostic_rhs[lev]->array(mfi); + out_flux = diagnostic_flux[lev]->array(mfi); + } amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) { - next(i,j,k) = eta(i,j,k) + dt * ConservativeRHS(eta, mu, i, j, k, dx.data(), L_local, floor); + const Set::Scalar rhs = ConservativeRHS(eta, mu, mask, i, j, k, dx.data(), L_local, floor); + const Set::Scalar raw = eta(i,j,k) + dt * rhs; + next(i,j,k) = clipping_local ? std::max(-1.0, std::min(1.0, raw)) : raw; + if (diagnostics_local) + { + const Set::Scalar alpha = std::max(0.0, std::min(1.0, 0.5*(eta(i,j,k)+1.0))); + const Set::Scalar center_m = mask(i,j,k)*(L_local*alpha*alpha*(1.0-alpha)*(1.0-alpha)+floor); + out_m(i,j,k) = center_m; out_mask(i,j,k) = mask(i,j,k); + out_raw(i,j,k) = raw; out_delta(i,j,k) = next(i,j,k)-eta(i,j,k); out_rhs(i,j,k) = rhs; + CellFlux(eta, mu, mask, i, j, k, dx.data(), L_local, floor, out_flux); + } }); } etanew_mf[lev]->FillBoundary(geom[lev].periodicity()); @@ -208,13 +312,40 @@ CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar using Operator::Spectral::FFT; FFT fft(geom[lev]); amrex::Vector hierarchy(lev + 1); - for (int ilev = 0; ilev <= lev; ++ilev) hierarchy[ilev] = etanew_mf[ilev].get(); + amrex::Vector> old_alpha(lev + 1); + for (int ilev = 0; ilev <= lev; ++ilev) + { + hierarchy[ilev] = etanew_mf[ilev].get(); + if (output_diagnostics) + { + old_alpha[ilev] = std::make_unique(etanew_mf[ilev]->boxArray(), etanew_mf[ilev]->DistributionMap(), 1, 0); + amrex::MultiFab::Copy(*old_alpha[ilev], *etanew_mf[ilev], 0, 0, 1, 0); + } + } // Form derivatives on a complete grid so coarse/fine patch edges have valid data. auto eta_mf = FFT::CompositeToUniform(hierarchy, geom, refRatio(), lev, time, *bc, 1, 1); eta_mf->FillBoundary(geom[lev].periodicity()); amrex::MultiFab mu_mf(eta_mf->boxArray(), eta_mf->DistributionMap(), 1, 1); ChemicalPotential(*eta_mf, mu_mf, geom[lev], gamma); mu_mf.FillBoundary(geom[lev].periodicity()); + amrex::Vector> mask_hierarchy(lev + 1); + amrex::Vector mask_pointers(lev + 1); + for (int ilev = 0; ilev <= lev; ++ilev) + { + mask_hierarchy[ilev] = std::make_unique( + etanew_mf[ilev]->boxArray(), etanew_mf[ilev]->DistributionMap(), 1, 1); + FillMobilityMask(ilev, *mask_hierarchy[ilev]); + mask_hierarchy[ilev]->FillBoundary(geom[ilev].periodicity()); + mask_pointers[ilev] = mask_hierarchy[ilev].get(); + } + auto mobility_mask = FFT::CompositeToUniform(mask_pointers, geom, refRatio(), lev, time, *bc, 1, 1); + mobility_mask->FillBoundary(geom[lev].periodicity()); + std::unique_ptr diagnostic_mobility_mf, diagnostic_flux_mf; + if (output_diagnostics) + { + diagnostic_mobility_mf = std::make_unique(eta_mf->boxArray(), eta_mf->DistributionMap(), 1, 0); + diagnostic_flux_mf = std::make_unique(eta_mf->boxArray(), eta_mf->DistributionMap(), AMREX_SPACEDIM, 0); + } amrex::Vector potential_hierarchy(lev + 1); for (int ilev = 0; ilev <= lev; ++ilev) potential_hierarchy[ilev] = intermediate[ilev].get(); FFT::UniformToComposite(mu_mf, potential_hierarchy, geom, refRatio(), lev, 1); @@ -222,17 +353,47 @@ CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar const auto dx = geom[lev].CellSizeArray(); const Set::Scalar L_local = L; const Set::Scalar floor = mobility_floor; + const bool diagnostics_local = output_diagnostics; for (amrex::MFIter mfi(*eta_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { const auto eta = eta_mf->const_array(mfi); const auto mu = mu_mf.const_array(mfi); + const auto mask = mobility_mask->const_array(mfi); const auto rhs = rhs_mf.array(mfi); + amrex::Array4 diagnostic_mobility, diagnostic_flux; + if (diagnostics_local) + { + diagnostic_mobility = diagnostic_mobility_mf->array(mfi); + diagnostic_flux = diagnostic_flux_mf->array(mfi); + } amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) { // Face gradients and their adjoint divergence dissipate the discrete energy. - rhs(i,j,k) = ConservativeRHS(eta, mu, i, j, k, dx.data(), L_local, floor); + rhs(i,j,k) = ConservativeRHS(eta, mu, mask, i, j, k, dx.data(), L_local, floor); + if (diagnostics_local) + { + const Set::Scalar alpha = std::max(0.0, std::min(1.0, 0.5*(eta(i,j,k)+1.0))); + const Set::Scalar center_m = mask(i,j,k)*(L_local*alpha*alpha*(1.0-alpha)*(1.0-alpha)+floor); + diagnostic_mobility(i,j,k) = center_m; + CellFlux(eta, mu, mask, i, j, k, dx.data(), L_local, floor, diagnostic_flux); + } }); } + if (output_diagnostics) + { + amrex::Vector dmask(lev+1), dmob(lev+1), dflux(lev+1), drhs(lev+1); + for (int ilev = 0; ilev <= lev; ++ilev) + { + dmask[ilev] = diagnostic_mobility_mask[ilev].get(); + dmob[ilev] = diagnostic_mobility[ilev].get(); + dflux[ilev] = diagnostic_flux[ilev].get(); + drhs[ilev] = diagnostic_rhs[ilev].get(); + } + FFT::UniformToComposite(*mobility_mask, dmask, geom, refRatio(), lev, 1); + FFT::UniformToComposite(*diagnostic_mobility_mf, dmob, geom, refRatio(), lev, 1); + FFT::UniformToComposite(*diagnostic_flux_mf, dflux, geom, refRatio(), lev, AMREX_SPACEDIM); + FFT::UniformToComposite(rhs_mf, drhs, geom, refRatio(), lev, 1); + } auto eta_hat_mf = fft.MakeSpectralFab(); auto rhs_hat_mf = fft.MakeSpectralFab(); fft.Forward(*eta_mf, eta_hat_mf); @@ -261,7 +422,27 @@ CahnHilliard::AdvanceDegenerateSpectral (int lev, Set::Scalar time, Set::Scalar }); } fft.Backward(eta_hat_mf, *eta_mf); + if (output_diagnostics) + { + amrex::Vector raw_hierarchy(lev+1); + for (int ilev = 0; ilev <= lev; ++ilev) raw_hierarchy[ilev] = diagnostic_raw[ilev].get(); + FFT::UniformToComposite(*eta_mf, raw_hierarchy, geom, refRatio(), lev, 1); + } FFT::UniformToComposite(*eta_mf, hierarchy, geom, refRatio(), lev, 1); + if (clipping) + for (int ilev = 0; ilev <= lev; ++ilev) + for (amrex::MFIter mfi(*etanew_mf[ilev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto eta = etanew_mf[ilev]->array(mfi); + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) + { eta(i,j,k) = std::max(-1.0, std::min(1.0, eta(i,j,k))); }); + } + if (output_diagnostics) + for (int ilev = 0; ilev <= lev; ++ilev) + { + amrex::MultiFab::Copy(*diagnostic_delta[ilev], *etanew_mf[ilev], 0, 0, 1, 0); + amrex::MultiFab::Subtract(*diagnostic_delta[ilev], *old_alpha[ilev], 0, 0, 1, 0); + } for (int ilev = 0; ilev <= lev; ++ilev) etanew_mf[ilev]->FillBoundary(geom[ilev].periodicity()); } #else @@ -349,6 +530,36 @@ CahnHilliard::Initialize (int lev) ic->Initialize(lev,etanew_mf); if (method == "realspace") ic->Initialize(lev,etaold_mf); + InitializeDiagnostics(lev); +} + +void +CahnHilliard::InitializeDiagnostics(int lev) +{ + if (!output_diagnostics) return; + amrex::MultiFab mask(etanew_mf[lev]->boxArray(), etanew_mf[lev]->DistributionMap(), 1, 0); + FillMobilityMask(lev, mask); + const bool sdch = mobility == "singly_degenerate"; + const Set::Scalar L_local = L, floor = mobility_floor; + for (amrex::MFIter mfi(*etanew_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto eta = etanew_mf[lev]->const_array(mfi); + const auto m = mask.const_array(mfi); + const auto out_m = diagnostic_mobility[lev]->array(mfi); + const auto out_mask = diagnostic_mobility_mask[lev]->array(mfi); + const auto out_raw = diagnostic_raw[lev]->array(mfi); + const auto out_delta = diagnostic_delta[lev]->array(mfi); + const auto out_rhs = diagnostic_rhs[lev]->array(mfi); + const auto out_flux = diagnostic_flux[lev]->array(mfi); + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + const Set::Scalar alpha = sdch ? std::max(0.0, std::min(1.0, 0.5*(eta(i,j,k)+1.0))) : 0.5; + out_mask(i,j,k) = m(i,j,k); + out_m(i,j,k) = sdch ? m(i,j,k)*(L_local*alpha*alpha*(1.0-alpha)*(1.0-alpha)+floor) : L_local; + out_raw(i,j,k) = eta(i,j,k); out_delta(i,j,k) = 0.0; out_rhs(i,j,k) = 0.0; + AMREX_D_TERM(out_flux(i,j,k,0)=0.0;, out_flux(i,j,k,1)=0.0;, out_flux(i,j,k,2)=0.0;); + }); + } } diff --git a/tests/SinglyDegenerateCahnHilliard/Readme.rst b/tests/SinglyDegenerateCahnHilliard/Readme.rst index 766d9d5f8a..81b9b39b54 100644 --- a/tests/SinglyDegenerateCahnHilliard/Readme.rst +++ b/tests/SinglyDegenerateCahnHilliard/Readme.rst @@ -9,21 +9,29 @@ Laplacian for the chemical potential. The spectral method uses FFTs to apply a discrete biharmonic stabilizer to the increment, including Nyquist modes. It does not guarantee energy decrease for arbitrary timesteps. Constant mobility remains the default. Periodic boundary conditions are required. -The realspace timestep must satisfy the explicit fourth-order stability constraint; -neither method clips the solution. +The realspace timestep must satisfy the explicit fourth-order stability constraint. +Neither method clips the solution by default; optional clipping can alter the +mean, so it is disabled in the conservation checks. -The regression sections check a linearized Fourier mode on anisotropic domains, -nonzero mobility floors, stationary phases, conservation, and energy dissipation -in two and three dimensions. Grid-scale checkerboards must decay, and a -strongly varying mobility case must dissipate energy at small timesteps. MPI -sections exercise patch/rank boundaries, and the spectral AMR section exercises -a partially refined grid with regridding. +The 12 regression sections keep a compact set of distinct checks. In both two +and three dimensions, realspace and spectral updates are checked against the +same linearized Fourier-mode solution, which also checks their agreement and +mass conservation. A spectral case checks a nonzero mobility floor; realspace +and spectral cases keep a pure phase stationary, and a further realspace case +checks an initial value outside the pure-phase interval. Spectral checkerboards +must decay in two and three dimensions, a strongly varying mobility case must +dissipate energy at a small timestep, and a regridded spectral AMR case checks +partial fine-level coverage. The 2D spectral Fourier-mode case also uses custom +input and field names and enables diagnostics, checking the namespaced fields. Realspace AMR does not reflux coarse/fine fluxes; use a uniform grid when exact mass conservation is required. Spectral AMR assembles a uniform finest-level grid, which requires memory proportional to that entire grid. -For example, run the two-dimensional serial sections with: +The checker reads both yt native plotfiles and Alamo HDF5 plotfiles. For example, +run the two-dimensional serial sections with: .. code-block:: bash scripts/runtests.py --dim=2 --serial --fft tests/SinglyDegenerateCahnHilliard + +For an HDF5-enabled build, add :code:`--hdf5` to select its matching executable. diff --git a/tests/SinglyDegenerateCahnHilliard/input b/tests/SinglyDegenerateCahnHilliard/input index f32de6f8bf..f037a9295e 100644 --- a/tests/SinglyDegenerateCahnHilliard/input +++ b/tests/SinglyDegenerateCahnHilliard/input @@ -1,282 +1,83 @@ -#@ [2d-realspace-amr] -#@ dim = 2 -#@ check-file = amr-real -#@ args = method=realspace timestep=0.000001 amr.max_level=1 amr.nsubsteps=1 -#@ args = refinement_threshold=0.015 amr.regrid_int=2 max_step=8 -#@ args = eta.ic.expression.region0="0.15+0.1*exp(-40*((x-0.5)^2+(y-0.75)^2))" -#@ -#@ [2d-spectral-linear-odd] -#@ dim = 2 -#@ fft = true -#@ check-file = linear-odd -#@ args = method=spectral amr.n_cell=31 31 31 amr.blocking_factor=1 max_step=1 -#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(30*pi*y/1.5)" -#@ -#@ [2d-spectral-unstabilized] -#@ dim = 2 -#@ fft = true -#@ check-file = nonlinear -#@ args = method=spectral spectral_stabilization=0 max_step=40 -#@ #@ [2d-realspace-linear] #@ dim = 2 -#@ check-file = linear -#@ args = method=realspace -#@ args = mobility_floor=0.0 +#@ check-file = linear-realspace +#@ args = method=realspace max_step=1 #@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" -#@ args = max_step=1 -#@ -#@ [2d-realspace-linear-floor] -#@ dim = 2 -#@ check-file = linear-floor -#@ args = method=realspace -#@ args = mobility_floor=0.0125 -#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" -#@ args = max_step=1 -#@ -#@ [2d-realspace-nonlinear] -#@ dim = 2 -#@ check-file = nonlinear -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = max_step=40 -#@ -#@ [2d-realspace-stationary] -#@ dim = 2 -#@ check-file = stationary -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="-1" -#@ args = max_step=4 -#@ -#@ [2d-realspace-stationary-outside] -#@ dim = 2 -#@ check-file = stationary-outside -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="1.2" -#@ args = max_step=4 #@ #@ [2d-spectral-linear] #@ dim = 2 -#@ check-file = linear #@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" -#@ args = max_step=1 +#@ check-file = linear-spectral-renamed +#@ args = method=spectral input_name=custom field_name=ch.eta old_field_name=ch.eta_old intermediate_name=ch.int diagnostics=true clipping=false max_step=1 +#@ args = custom.ic.type=expression custom.ic.expression.constant.pi=3.14159265358979323846 +#@ args = custom.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" +#@ args = custom.bc.constant.type.xlo=periodic custom.bc.constant.type.xhi=periodic +#@ args = custom.bc.constant.type.ylo=periodic custom.bc.constant.type.yhi=periodic #@ #@ [2d-spectral-linear-floor] #@ dim = 2 -#@ check-file = linear-floor #@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0125 +#@ check-file = linear-spectral-floor +#@ args = method=spectral mobility_floor=0.0125 max_step=1 #@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)" -#@ args = max_step=1 -#@ -#@ [2d-spectral-nonlinear] -#@ dim = 2 -#@ check-file = nonlinear -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = max_step=40 -#@ -#@ [2d-spectral-stationary] -#@ dim = 2 -#@ check-file = stationary -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="-1" -#@ args = max_step=4 -#@ -#@ [2d-spectral-stationary-outside] -#@ dim = 2 -#@ check-file = stationary-outside -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="1.2" -#@ args = max_step=4 #@ #@ [3d-realspace-linear] #@ dim = 3 -#@ check-file = linear -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" -#@ args = max_step=1 -#@ -#@ [3d-realspace-linear-floor] -#@ dim = 3 -#@ check-file = linear-floor -#@ args = method=realspace -#@ args = mobility_floor=0.0125 +#@ check-file = linear-realspace +#@ args = method=realspace max_step=1 #@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" -#@ args = max_step=1 -#@ -#@ [3d-realspace-nonlinear] -#@ dim = 3 -#@ check-file = nonlinear -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = max_step=40 -#@ -#@ [3d-realspace-stationary] -#@ dim = 3 -#@ check-file = stationary -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="-1" -#@ args = max_step=4 -#@ -#@ [3d-realspace-stationary-outside] -#@ dim = 3 -#@ check-file = stationary-outside -#@ args = method=realspace -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="1.2" -#@ args = max_step=4 #@ #@ [3d-spectral-linear] #@ dim = 3 -#@ check-file = linear #@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 +#@ check-file = linear-spectral +#@ args = method=spectral max_step=1 #@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" -#@ args = max_step=1 -#@ -#@ [3d-spectral-linear-floor] -#@ dim = 3 -#@ check-file = linear-floor -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0125 -#@ args = eta.ic.expression.region0="0.15+1e-6*cos(2*pi*x)*cos(2*pi*y/1.5)*cos(2*pi*z/2)" -#@ args = max_step=1 -#@ -#@ [3d-spectral-nonlinear] -#@ dim = 3 -#@ check-file = nonlinear -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = max_step=40 -#@ -#@ [3d-spectral-stationary] -#@ dim = 3 -#@ check-file = stationary -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="-1" -#@ args = max_step=4 #@ -#@ [3d-spectral-stationary-outside] -#@ dim = 3 -#@ check-file = stationary-outside -#@ fft = true -#@ args = method=spectral -#@ args = mobility_floor=0.0 -#@ args = eta.ic.expression.region0="1.2" -#@ args = max_step=4 -#@ -#@ [2d-realspace-mpi] +#@ [2d-realspace-pure-phase] #@ dim = 2 -#@ nprocs = 2 -#@ check-file = nonlinear -#@ args = method=realspace -#@ args = max_step=40 +#@ check-file = stationary +#@ args = method=realspace mobility_floor=0 eta.ic.expression.region0="-1" max_step=4 #@ -#@ [2d-spectral-mpi] +#@ [2d-spectral-pure-phase] #@ dim = 2 -#@ nprocs = 2 -#@ check-file = nonlinear -#@ args = method=spectral -#@ args = max_step=40 #@ fft = true +#@ check-file = stationary +#@ args = method=spectral mobility_floor=0 eta.ic.expression.region0="-1" max_step=4 #@ -#@ [2d-spectral-amr] +#@ [2d-realspace-outside-phase] #@ dim = 2 -#@ fft = true -#@ check-file = amr -#@ args = method=spectral amr.max_level=1 amr.nsubsteps=1 -#@ args = refinement_threshold=0.015 amr.regrid_int=2 max_step=8 -#@ args = eta.ic.expression.region0="0.15+0.1*exp(-40*((x-0.5)^2+(y-0.75)^2))" - -#@ [2d-spectral-checkerboard-x] +#@ check-file = stationary-outside +#@ args = method=realspace mobility_floor=0 eta.ic.expression.region0="1.2" max_step=4 +#@ +#@ [2d-spectral-checkerboard] #@ dim = 2 #@ fft = true #@ check-file = checkerboard-x #@ args = method=spectral max_step=4 #@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)" #@ -#@ [2d-spectral-checkerboard-y] -#@ dim = 2 -#@ fft = true -#@ check-file = checkerboard-y -#@ args = method=spectral max_step=4 -#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*y/1.5)" -#@ -#@ [2d-spectral-checkerboard-xy] -#@ dim = 2 -#@ fft = true -#@ check-file = checkerboard-xy -#@ args = method=spectral max_step=4 -#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)*sin(32*pi*y/1.5)" -#@ -#@ [2d-spectral-strong-mobility] -#@ dim = 2 -#@ args = timestep=1e-8 -#@ -#@ [2d-spectral-strong-mobility-small] -#@ dim = 2 -#@ args = timestep=1e-10 -#@ -#@ [3d-spectral-checkerboard-x] +#@ [3d-spectral-checkerboard] #@ dim = 3 #@ fft = true #@ check-file = checkerboard-x #@ args = method=spectral max_step=4 #@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)" #@ -#@ [3d-spectral-checkerboard-y] -#@ dim = 3 -#@ fft = true -#@ check-file = checkerboard-y -#@ args = method=spectral max_step=4 -#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*y/1.5)" -#@ -#@ [3d-spectral-checkerboard-z] -#@ dim = 3 -#@ fft = true -#@ check-file = checkerboard-z -#@ args = method=spectral max_step=4 -#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*z/2)" -#@ -#@ [3d-spectral-checkerboard-xyz] -#@ dim = 3 -#@ fft = true -#@ check-file = checkerboard-xyz -#@ args = method=spectral max_step=4 -#@ args = eta.ic.expression.region0="0.15+0.01*sin(32*pi*x)*sin(32*pi*y/1.5)*sin(32*pi*z/2)" -#@ -#@ [3d-spectral-strong-mobility] -#@ dim = 3 -#@ args = timestep=1e-8 -#@ -#@ [3d-spectral-strong-mobility-small] -#@ dim = 3 -#@ args = timestep=1e-10 -#@ -#@ [*spectral-strong-mobility*] +#@ [2d-spectral-strong-mobility] +#@ dim = 2 #@ fft = true #@ check-file = strong-mobility -#@ args = method=spectral amr.n_cell=16 16 16 max_step=1 +#@ args = method=spectral timestep=1e-8 amr.n_cell=16 16 max_step=1 #@ args = eta.ic.expression.region0="(0.17623171536106855)*(x>=0 and x<0.0625)+(0.9999915189178985)*(x>=0.0625 and x<0.125)+(0.99957222732831352)*(x>=0.125 and x<0.1875)+(-0.33828120309429632)*(x>=0.1875 and x<0.25)+(0.96308521781430201)*(x>=0.25 and x<0.3125)+(0.86409632475355058)*(x>=0.3125 and x<0.375)+(0.98990644052604582)*(x>=0.375 and x<0.4375)+(-0.98167620233468933)*(x>=0.4375 and x<0.5)+(0.99984268465982185)*(x>=0.5 and x<0.5625)+(0.99898999623396001)*(x>=0.5625 and x<0.625)+(-0.54394735059310761)*(x>=0.625 and x<0.6875)+(-0.99762312075226711)*(x>=0.6875 and x<0.75)+(-0.99954221397642695)*(x>=0.75 and x<0.8125)+(0.99930823691403747)*(x>=0.8125 and x<0.875)+(0.99991743650726583)*(x>=0.875 and x<0.9375)+(0.89293086073486305)*(x>=0.9375 and x<1)" +#@ +#@ [2d-spectral-amr] +#@ dim = 2 +#@ fft = true +#@ check-file = amr +#@ args = method=spectral amr.max_level=1 amr.nsubsteps=1 +#@ args = refinement_threshold=0.015 amr.regrid_int=2 max_step=8 +#@ args = eta.ic.expression.region0="0.15+0.1*exp(-40*((x-0.5)^2+(y-0.75)^2))" alamo.program = cahnhilliard plot_file = tests/SinglyDegenerateCahnHilliard/output diff --git a/tests/SinglyDegenerateCahnHilliard/test b/tests/SinglyDegenerateCahnHilliard/test index b12770a396..958ad624b3 100755 --- a/tests/SinglyDegenerateCahnHilliard/test +++ b/tests/SinglyDegenerateCahnHilliard/test @@ -1,7 +1,9 @@ #!/usr/bin/env python3 """Check conservation, dissipation, pure phases, and Fourier-mode evolution.""" +import re import sys from pathlib import Path +from types import SimpleNamespace import numpy as np import yt @@ -9,21 +11,122 @@ import yt yt.set_log_level(50) outdir = Path(sys.argv[1]) case = sys.argv[2] -plots = sorted(outdir.glob("*cell")) +plots = list(outdir.glob("*cell")) + list(outdir.glob("*cell.h5")) +plots.sort(key=lambda path: int(re.search(r"(\d+)cell", path.name).group(1))) assert len(plots) >= 2, "Missing initial/final plotfiles" +eta_field = "ch.eta" if case == "linear-spectral-renamed" else "eta" +int_field = "ch.int" if case == "linear-spectral-renamed" else "int" + + def read(path): + if path.suffix == ".h5": + return read_hdf5(path, eta_field) ds = yt.load(str(path)) level = ds.index.max_level dimensions = ds.domain_dimensions.copy() dimensions[:ds.dimensionality] *= 2**level grid = ds.covering_grid(level, ds.domain_left_edge, dimensions) - eta = np.asarray(grid[("boxlib", "eta")]).squeeze() + eta = np.asarray(grid[("boxlib", eta_field)]).squeeze() dx = np.asarray(ds.domain_width)[:eta.ndim] / np.array(eta.shape) return ds, eta, dx +def read_hdf5(path, eta_name): + """Read Alamo's AMReX cell plotfiles without requiring yt's native reader.""" + import h5py + + with h5py.File(path, "r") as handle: + components = {} + for key, value in handle.attrs.items(): + if key.startswith("component_"): + name = value.decode() if isinstance(value, bytes) else str(value) + components[name] = int(key.rsplit("_", 1)[1]) + level_names = sorted( + (name for name in handle if name.startswith("level_")), + key=lambda name: int(name.split("_")[-1]), + ) + max_level = int(level_names[-1].split("_")[-1]) + def metadata(name, default=None): + if name in handle.attrs: + return handle.attrs[name] + if name in handle["level_0"].attrs: + return handle["level_0"].attrs[name] + if name in handle: + return np.asarray(handle[name]) + if default is not None: + return default + raise KeyError(f"Missing HDF5 plotfile metadata: {name}") + + ndim = int(np.asarray(metadata("dim")).reshape(-1)[0]) + dx = np.asarray(metadata("Vec_dx"), dtype=float)[:ndim] + prob_lo = np.asarray(metadata("prob_lo", np.zeros(ndim)), dtype=float)[:ndim] + time = float(np.asarray(metadata("time")).reshape(-1)[0]) + + def box_bounds(box): + names = box.dtype.names + if names: + axes = "ijk"[:ndim] + lo = np.array([box[f"lo_{axis}"] for axis in axes], dtype=int) + hi = np.array([box[f"hi_{axis}"] for axis in axes], dtype=int) + return lo, hi + values = np.asarray(box, dtype=int) + return values[:ndim], values[ndim:] + + first_boxes = handle["level_0"]["boxes"][:] + domain_dimensions = np.max([box_bounds(box)[1] for box in first_boxes], axis=0) + 1 + finest_shape = tuple((domain_dimensions * 2**max_level).tolist()) + eta_component = next((index for name, index in components.items() + if name == eta_name or name.endswith("." + eta_name.split(".")[-1])), None) + if eta_component is None: + raise KeyError(f"HDF5 plotfile has no {eta_name!r} component; found {list(components)}") + eta_key = next(name for name, index in components.items() if index == eta_component) + fields = {name: np.full(finest_shape, np.nan) for name in components} + grids = [] + for level_name in level_names: + level = int(level_name.split("_")[-1]) + group = handle[level_name] + boxes = group["boxes"][:] + offsets = np.asarray(group["data:offsets=0"]) + values = np.asarray(group["data:datatype=0"]) + bounds = [box_bounds(box) for box in boxes] + grids.extend(SimpleNamespace(Level=level, ActiveDimensions=hi - lo + 1) + for lo, hi in bounds) + for box_index, box in enumerate(boxes): + lo, hi = bounds[box_index] + box_shape = hi - lo + 1 + cell_count = int(np.prod(box_shape)) + begin = int(offsets[box_index]) + scale = 2**(max_level - level) + starts, stops = lo * scale, (hi + 1) * scale + # Coarse AMR data fills its represented fine cells; finer levels overwrite it. + for name, component in components.items(): + raw = values[begin + component * cell_count: + begin + (component + 1) * cell_count] + patch = raw.reshape(box_shape[::-1]).transpose(tuple(range(ndim - 1, -1, -1))) + for axis in range(ndim): + patch = np.repeat(patch, scale, axis=axis) + fields[name][tuple(slice(a, b) for a, b in zip(starts, stops))] = patch + + width = dx * domain_dimensions + ds = SimpleNamespace( + current_time=time, + domain_dimensions=domain_dimensions, + domain_width=width, + domain_left_edge=prob_lo, + dimensionality=ndim, + index=SimpleNamespace(max_level=max_level, grids=grids), + mean_value=float(np.nanmean(fields[eta_key])), + fields=fields, + covering_grid=lambda level, *_args, **_kwargs: { + ("boxlib", name): field[tuple(slice(None, None, 2**(max_level-level)) for _ in range(ndim))] + for name, field in fields.items() + }, + ) + return ds, fields[eta_key], dx / 2**max_level + + def energy(eta, dx): density = 0.25 * (eta**2 - 1)**2 for axis, h in enumerate(dx): @@ -32,16 +135,39 @@ def energy(eta, dx): def mean(ds): + if hasattr(ds, "mean_value"): + return ds.mean_value cells = ds.all_data() volumes = cells[("index", "cell_volume")] - return float((cells[("boxlib", "eta")] * volumes).sum() / volumes.sum()) + return float((cells[("boxlib", eta_field)] * volumes).sum() / volumes.sum()) initial_ds, initial, dx = read(plots[0]) final_ds, final, _ = read(plots[-1]) assert np.isfinite(final).all(), "Nonfinite SDCH solution" -if case != "amr-real": - assert abs(mean(final_ds) - mean(initial_ds)) < 2e-12, "SDCH changed the mean" +if case == "linear-spectral-renamed": + diagnostic_names = ("ch.mobility", "ch.rhs", "ch.alpha_raw", "ch.alpha_delta", + "ch.mobility_mask", "ch.flux_x", "ch.flux_y") + if hasattr(final_ds, "fields"): + diagnostic_fields = final_ds.fields + else: + grid = final_ds.covering_grid(final_ds.index.max_level, final_ds.domain_left_edge, + final_ds.domain_dimensions * 2**final_ds.index.max_level) + diagnostic_fields = {name: np.asarray(grid[("boxlib", name)]).squeeze() + for name in diagnostic_names} + for name in diagnostic_names: + assert name in diagnostic_fields, f"Missing diagnostic field {name}" + assert np.isfinite(diagnostic_fields[name]).all(), f"Nonfinite diagnostic {name}" + alpha = np.clip(0.5 * (initial + 1.0), 0.0, 1.0) + expected_mobility = alpha**2 * (1.0 - alpha)**2 + np.testing.assert_allclose(diagnostic_fields["ch.mobility"], expected_mobility, + rtol=0, atol=3e-16) + np.testing.assert_array_equal(diagnostic_fields["ch.mobility_mask"], 1.0) + np.testing.assert_allclose(diagnostic_fields["ch.alpha_raw"], final, rtol=0, atol=2e-14) + np.testing.assert_allclose(diagnostic_fields["ch.alpha_delta"], final - initial, + rtol=0, atol=2e-14) + assert np.max(np.abs(diagnostic_fields["ch.rhs"])) > 1e-10, "Diagnostic RHS is zero" +assert abs(mean(final_ds) - mean(initial_ds)) < 2e-12, "SDCH changed the mean" assert float(final_ds.current_time) > float(initial_ds.current_time) if case.startswith("stationary"): @@ -54,16 +180,14 @@ elif case.startswith("linear"): mobility = alpha**2 * (1-alpha)**2 + floor widths = np.asarray(initial_ds.domain_width)[:initial.ndim] modes = np.ones(initial.ndim) - if case == "linear-odd": - modes[1] = 15 wave = 2*np.pi*modes/widths lap = np.sum(4*np.sin(wave*dx/2)**2/dx**2) curvature = 3*base**2 - 1 + 0.01*lap potential = np.asarray(final_ds.covering_grid(0, final_ds.domain_left_edge, - final_ds.domain_dimensions)[("boxlib", "int")]).squeeze() + final_ds.domain_dimensions)[("boxlib", int_field)]).squeeze() expected_potential = base**3 - base + (initial-base) * curvature np.testing.assert_allclose(potential, expected_potential, rtol=0, atol=1e-12) - spectral = "spectral" in outdir.name + spectral = "spectral" in case rate = -mobility * lap * curvature if spectral: stabilizer = 0.01 * (1/16 + floor)