#include <ERF_NodalReconstruction.H>

Collaboration diagram for NodalReconstruction:

[legend]

Public Types
using	Real = amrex::Real

Public Member Functions
	NodalReconstruction (amrex::Box const &face_box, amrex::Geometry const &geom)

amrex::Box const &	nodalBox () const noexcept

amrex::FArrayBox	makeReference (amrex::FArrayBox const &T) const

std::pair< amrex::FArrayBox, SolveInfo >	solve (amrex::FArrayBox const &T, amrex::FArrayBox const &reference, VariationOperator var_op, Real relative_tolerance=Real(1.e-10), Real smoothing=Real(1.e-2), Real min_regularization=Real(2.5e-4), int max_iterations=1000)

Real	maxAverageError (amrex::FArrayBox const &T, amrex::FArrayBox const &S) const
	max_ij \| avg4(S) - T \| More...

Real	totalSquaredVariation (amrex::FArrayBox const &S) const

Private Member Functions
void	applyA (amrex::FArrayBox const &S, amrex::FArrayBox &T) const
	(A S)(i,j) = 1/4 [ S(i,j) + S(i+1,j) + S(i,j+1) + S(i+1,j+1) ] More...

void	applyAT (amrex::FArrayBox const &T, amrex::FArrayBox &R) const
	Adjoint of applyA: scatter each cell value to its four nodes. More...

void	applyD (amrex::FArrayBox const &S, amrex::FArrayBox &difx, amrex::FArrayBox &dify) const

void	applyDT (amrex::FArrayBox const &difx, amrex::FArrayBox const &dify, amrex::FArrayBox &R) const

void	applyDTD (amrex::FArrayBox const &S, amrex::FArrayBox &R) const

void	applyL (amrex::FArrayBox const &S, amrex::FArrayBox &lap) const

void	applyLT (amrex::FArrayBox const &lap, amrex::FArrayBox &R) const

void	applyLTL (amrex::FArrayBox const &S, amrex::FArrayBox &R) const

void	applyVariationOperator (amrex::FArrayBox const &S, amrex::FArrayBox &R) const

void	applyM (amrex::FArrayBox const &S, amrex::FArrayBox &R) const
	R = ( A^T A + lambda V^T V + mu I ) S. More...

void	computeDiagonal (amrex::FArrayBox &diag) const

void	computeDiagnostics (amrex::FArrayBox const &T, amrex::FArrayBox const &S, SolveInfo &info) const

Real	maxAbs (amrex::FArrayBox const &F) const

SolveInfo	conjugateGradient (amrex::FArrayBox const &rhs, amrex::FArrayBox &x, Real relative_tolerance, int max_iterations) const

void	applyJacobi (amrex::FArrayBox const &diag, amrex::FArrayBox const &r, amrex::FArrayBox &z) const

Real	dot (amrex::FArrayBox const &a, amrex::FArrayBox const &b) const

Private Attributes
amrex::Box	m_face_box

amrex::Box	m_nodal_box

int	m_ilo =0

int	m_ihi =-1

int	m_jlo =0

int	m_jhi =-1

amrex::Real	m_gx = one

amrex::Real	m_gy = one

amrex::Real	m_gxx = one

amrex::Real	m_gyy = one

amrex::Real	m_lambda = zero

amrex::Real	m_mu = zero

VariationOperator	m_var_op = VariationOperator::FirstDeriv

amrex::FArrayBox	m_cc_scratch

amrex::FArrayBox	m_nd_scratch

Static Private Attributes
static constexpr amrex::Real	deviation_factor = amrex::Real(2.0)

static constexpr amrex::Real	relief_factor = amrex::Real(0.25)
	... nor than this fraction of the relief of the data. More...

static constexpr amrex::Real	mu_growth = amrex::Real(4.0)

static constexpr int	max_attempts = 24

Detailed Description

Reconstruct a nodal field S from data T that lives at cell centers in (x,y), such that the four-node average of S matches T as closely as a bounded, smooth nodal field can.

Why we do not simply invert the averaging operator

The four-node averaging operator A,

(A S)(i,j) = 1/4 [ S(i,j) + S(i+1,j) + S(i,j+1) + S(i+1,j+1) ] ,

has symbol cos(kx/2) cos(ky/2), which vanishes at the grid Nyquist mode. Inverting A exactly – e.g. with the back substitution

S(i+1,j+1) = 4 T(i,j) - S(i,j) - S(i+1,j) - S(i,j+1)

– therefore amplifies grid-scale content of T without bound. The Green's function of that recursion is 4 / [(1+x)(1+y)], a doubly alternating cumulative sum, so a checkerboard component of T of amplitude a produces a nodal response of amplitude a*i*j. On a 128^2 grid an O(10 m) grid-scale roughness in the WRF terrain thus yields nodal heights that are kilometers away from the WRF terrain, including physically impossible negative elevations (ERF issue #3709). Adding a free homogeneous solution cannot repair this: the null space of A is spanned by (-1)^i f(j) + (-1)^j g(i), which cannot represent the (-1)^(i+j) i j growth of the particular solution. The failure is silent, because checking that avg4(S) reproduces T is then satisfied by construction.

What we do instead

Let S_ref be the direct interpolation of T to the nodes (makeReference) and R = T - A S_ref the part of T that interpolation fails to reproduce. We solve for a correction E to the interpolant,

min_E  || A E - R ||^2 + lambda || V E ||^2 + mu || E ||^2 ,   S = S_ref + E,

where V is a roughness operator (first derivative or Laplacian). The normal equations

( A^T A + lambda V^T V + mu I ) E = A^T R

are symmetric positive definite for any mu > 0, so E exists, is unique, and – unlike the exact inverse – is bounded. Specifically, writing ( A^T A + lambda V^T V + mu I ) E = A^T R and pairing with E gives

|| A E ||^2 + lambda || V E ||^2 + mu || E ||^2 <= || A E || || R || ,

and since 2 sqrt(mu) ||A E|| ||E|| <= ||A E||^2 + mu ||E||^2, we get the a priori bound

|| E || <= || R || / ( 2 sqrt(mu) ) .

The correction is therefore controlled by how badly interpolation already fails, never by the conditioning of A.

Choosing mu

A single fixed mu cannot serve both smooth and rough terrain: small mu recovers well-resolved features almost exactly but lets rough fields drift far from the interpolant, and large mu does the reverse. We therefore start from a small mu and raise it until the correction satisfies the pointwise cap

max |E| <= min( 2 max|R| , 1/4 (max T - min T) ) ,

i.e. the nodes may not move away from the direct interpolation by more than twice the interpolation misfit, nor by more than a quarter of the relief of the data. Because ||E|| -> 0 as mu -> infinity the loop always terminates, and on failure we fall back to the interpolant itself. Smooth terrain keeps the smallest mu and is reproduced to well below the interpolation error; rough terrain is pulled back toward the interpolant instead of blowing up.

The averaging error max |avg4(S) - T| returned in SolveInfo is now a genuine diagnostic – it is no longer satisfied by construction – and is reported alongside the same quantity for the plain interpolant and the range of S, so callers can sanity check the result.

Member Typedef Documentation

◆ Real

using NodalReconstruction::Real = amrex::Real

Constructor & Destructor Documentation

◆ NodalReconstruction()

NodalReconstruction::NodalReconstruction	(	amrex::Box const &	face_box,
		amrex::Geometry const &	geom
	)

inlineexplicit

         : m_face_box(face_box),
           m_ilo(face_box.smallEnd(0)),
           m_ihi(face_box.bigEnd(0)),
           m_jlo(face_box.smallEnd(1)),
           m_jhi(face_box.bigEnd(1))
     {
         if (face_box.length(2) != 1) {
             throw std::invalid_argument(
                 "face_box must contain exactly one z cell.");
         }
         if (face_box.smallEnd(2) != 0) {
             throw std::invalid_argument(
                 "face_box must be a slab at k = 0; all fabs are indexed at k = 0.");
         }
  
         m_nodal_box = amrex::surroundingNodes(face_box);
         m_nodal_box.setSmall(2, 0);
         m_nodal_box.setBig(2, 0);
  
         Real const dx = geom.CellSizeArray()[0];
         Real const dy = geom.CellSizeArray()[1];
         if (!(dx > zero) || !(dy > zero)) {
             throw std::invalid_argument("Cell sizes must be positive.");
         }
  
         // Non-dimensional, O(1) stencil weights that respect grid anisotropy.
         // They are unity on an isotropic grid, so the smoothing weight has the
         // same meaning regardless of the mesh spacing.
         Real const h = std::min(dx,dy);
         m_gx  = h / dx;
         m_gy  = h / dy;
         m_gxx = m_gx * m_gx;
         m_gyy = m_gy * m_gy;
     }

Member Function Documentation

◆ applyA()

void NodalReconstruction::applyA	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	T
	)		const

inlineprivate

(A S)(i,j) = 1/4 [ S(i,j) + S(i+1,j) + S(i,j+1) + S(i+1,j+1) ]

     {
         auto const S_arr = S.const_array();
         auto const T_arr = T.array();
         for (int j=m_jlo; j<=m_jhi; ++j) {
             for (int i=m_ilo; i<=m_ihi; ++i) {
                 T_arr(i,j,0) = fourth * ( S_arr(i  ,j  ,0) + S_arr(i+1,j  ,0)
                                         + S_arr(i  ,j+1,0) + S_arr(i+1,j+1,0) );
             }
         }
     }

Referenced by applyM(), and solve().

Here is the caller graph for this function:

◆ applyAT()

void NodalReconstruction::applyAT	(	amrex::FArrayBox const &	T,
		amrex::FArrayBox &	R
	)		const

inlineprivate

Adjoint of applyA: scatter each cell value to its four nodes.

     {
         R.setVal<amrex::RunOn::Host>(zero);
         auto const T_arr = T.const_array();
         auto const R_arr = R.array();
         for (int j=m_jlo; j<=m_jhi; ++j) {
             for (int i=m_ilo; i<=m_ihi; ++i) {
                 Real const q = fourth * T_arr(i,j,0);
                 R_arr(i  ,j  ,0) += q;
                 R_arr(i+1,j  ,0) += q;
                 R_arr(i  ,j+1,0) += q;
                 R_arr(i+1,j+1,0) += q;
             }
         }
     }

Referenced by applyM(), and solve().

Here is the caller graph for this function:

◆ applyD()

void NodalReconstruction::applyD	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	difx,
		amrex::FArrayBox &	dify
	)		const

inlineprivate

     {
         difx.setVal<amrex::RunOn::Host>(zero);
         dify.setVal<amrex::RunOn::Host>(zero);
  
         auto const S_arr    = S.const_array();
         auto const difx_arr = difx.array();
         auto const dify_arr = dify.array();
  
         for (int j = m_jlo; j <= m_jhi + 1; ++j) {
             for (int i = m_ilo; i <= m_ihi; ++i) {
                 difx_arr(i, j, 0) = m_gx * ( S_arr(i+1, j, 0) - S_arr(i, j, 0) );
             }
         }
  
         for (int j = m_jlo; j <= m_jhi; ++j) {
             for (int i = m_ilo; i <= m_ihi + 1; ++i) {
                 dify_arr(i, j, 0) = m_gy * ( S_arr(i, j+1, 0) - S_arr(i, j, 0) );
             }
         }
     }

Referenced by applyDTD().

Here is the caller graph for this function:

◆ applyDT()

void NodalReconstruction::applyDT	(	amrex::FArrayBox const &	difx,
		amrex::FArrayBox const &	dify,
		amrex::FArrayBox &	R
	)		const

inlineprivate

     {
         R.setVal<amrex::RunOn::Host>(zero);
  
         auto const qx_arr = difx.const_array();
         auto const qy_arr = dify.const_array();
         auto const R_arr  = R.array();
  
         for (int j = m_jlo; j <= m_jhi + 1; ++j) {
             for (int i = m_ilo; i <= m_ihi; ++i) {
                 R_arr(i  , j, 0) -= m_gx * qx_arr(i, j, 0);
                 R_arr(i+1, j, 0) += m_gx * qx_arr(i, j, 0);
             }
         }
  
         for (int j = m_jlo; j <= m_jhi; ++j) {
             for (int i = m_ilo; i <= m_ihi + 1; ++i) {
                 R_arr(i, j  , 0) -= m_gy * qy_arr(i, j, 0);
                 R_arr(i, j+1, 0) += m_gy * qy_arr(i, j, 0);
             }
         }
     }

Referenced by applyDTD().

Here is the caller graph for this function:

◆ applyDTD()

void NodalReconstruction::applyDTD	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	R
	)		const

inlineprivate

     {
         amrex::Box x_edge_box = m_nodal_box;
         amrex::Box y_edge_box = m_nodal_box;
         x_edge_box.growHi(0,-1);
         y_edge_box.growHi(1,-1);
         amrex::FArrayBox difx(x_edge_box,1,amrex::The_Managed_Arena());
         amrex::FArrayBox dify(y_edge_box,1,amrex::The_Managed_Arena());
         applyD(S, difx, dify);
         applyDT(difx, dify, R);
     }

Referenced by applyVariationOperator().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ applyJacobi()

void NodalReconstruction::applyJacobi	(	amrex::FArrayBox const &	diag,
		amrex::FArrayBox const &	r,
		amrex::FArrayBox &	z
	)		const

inlineprivate

     {
         auto const d_arr = diag.const_array();
         auto const r_arr = r.const_array();
         auto const z_arr = z.array();
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 z_arr(i,j,0) = r_arr(i,j,0) / d_arr(i,j,0);
             }
         }
     }

Referenced by conjugateGradient().

Here is the caller graph for this function:

◆ applyL()

void NodalReconstruction::applyL	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	lap
	)		const

inlineprivate

     {
         lap.setVal<amrex::RunOn::Host>(zero);
  
         auto const S_arr   = S.const_array();
         auto const lap_arr = lap.array();
  
         for (int j = m_jlo + 1; j <= m_jhi; ++j) {
             for (int i = m_ilo + 1; i <= m_ihi; ++i) {
                 lap_arr(i, j, 0) = m_gxx * (         S_arr(i+1, j  , 0)
                                             - two *  S_arr(i  , j  , 0)
                                                    + S_arr(i-1, j  , 0) )
                                  + m_gyy * (         S_arr(i  , j+1, 0)
                                             - two *  S_arr(i  , j  , 0)
                                                    + S_arr(i  , j-1, 0) );
             }
         }
     }

Referenced by applyLTL().

Here is the caller graph for this function:

◆ applyLT()

void NodalReconstruction::applyLT	(	amrex::FArrayBox const &	lap,
		amrex::FArrayBox &	R
	)		const

inlineprivate

     {
         R.setVal<amrex::RunOn::Host>(zero);
  
         auto const q_arr = lap.const_array();
         auto const R_arr = R.array();
  
         for (int j = m_jlo + 1; j <= m_jhi; ++j) {
             for (int i = m_ilo + 1; i <= m_ihi; ++i) {
                 const amrex::Real val = q_arr(i, j, 0);
                 R_arr(i-1, j  , 0) += m_gxx * val;
                 R_arr(i+1, j  , 0) += m_gxx * val;
                 R_arr(i  , j-1, 0) += m_gyy * val;
                 R_arr(i  , j+1, 0) += m_gyy * val;
                 R_arr(i  , j  , 0) -= two * (m_gxx + m_gyy) * val;
             }
         }
     }

Referenced by applyLTL().

Here is the caller graph for this function:

◆ applyLTL()

void NodalReconstruction::applyLTL	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	R
	)		const

inlineprivate

     {
         amrex::Box interior_nodal_box = m_nodal_box;
         interior_nodal_box.grow(0,-1);
         interior_nodal_box.grow(1,-1);
         if (!interior_nodal_box.ok()) {
             // Fewer than two cells in a direction: there is no interior node to
             // evaluate the Laplacian at.  mu keeps the operator definite.
             R.setVal<amrex::RunOn::Host>(zero);
             return;
         }
         amrex::FArrayBox lap(interior_nodal_box,1,amrex::The_Managed_Arena());
         applyL(S, lap);
         applyLT(lap, R);
     }

Referenced by applyVariationOperator().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ applyM()

void NodalReconstruction::applyM	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	R
	)		const

inlineprivate

R = ( A^T A + lambda V^T V + mu I ) S.

     {
         applyA(S,m_cc_scratch);
         applyAT(m_cc_scratch,R);
         applyVariationOperator(S,m_nd_scratch);
  
         auto const R_arr = R.array();
         auto const v_arr = m_nd_scratch.const_array();
         auto const S_arr = S.const_array();
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 R_arr(i,j,0) += m_lambda * v_arr(i,j,0) + m_mu * S_arr(i,j,0);
             }
         }
     }

Referenced by conjugateGradient().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ applyVariationOperator()

void NodalReconstruction::applyVariationOperator	(	amrex::FArrayBox const &	S,
		amrex::FArrayBox &	R
	)		const

inlineprivate

     {
         switch (m_var_op) {
         case VariationOperator::FirstDeriv:
             applyDTD(S, R);
             break;
  
         case VariationOperator::Laplacian:
             applyLTL(S, R);
             break;
  
         default:
             throw std::runtime_error("Unsupported variation operator.");
         }
     }

Referenced by applyM().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ computeDiagnostics()

void NodalReconstruction::computeDiagnostics	(	amrex::FArrayBox const &	T,
		amrex::FArrayBox const &	S,
		SolveInfo &	info
	)		const

inlineprivate

     {
         auto const S_arr = S.const_array();
         Real smin =  std::numeric_limits<Real>::max();
         Real smax = -std::numeric_limits<Real>::max();
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 smin = std::min(smin, S_arr(i,j,0));
                 smax = std::max(smax, S_arr(i,j,0));
             }
         }
         info.min_value = smin;
         info.max_value = smax;
         info.max_average_error = maxAverageError(T,S);
     }

Referenced by solve().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ computeDiagonal()

void NodalReconstruction::computeDiagonal ( amrex::FArrayBox & diag ) const

inlineprivate

Diagonal of M, assembled by accumulating the square of every stencil coefficient that touches a node. Used as a Jacobi preconditioner.

     {
         diag.setVal<amrex::RunOn::Host>(m_mu);
         auto const d_arr = diag.array();
  
         // A^T A: each incident cell contributes (1/4)^2
         Real const qa = fourth * fourth;
         for (int j=m_jlo; j<=m_jhi; ++j) {
             for (int i=m_ilo; i<=m_ihi; ++i) {
                 d_arr(i  ,j  ,0) += qa;
                 d_arr(i+1,j  ,0) += qa;
                 d_arr(i  ,j+1,0) += qa;
                 d_arr(i+1,j+1,0) += qa;
             }
         }
  
         // lambda V^T V
         if (m_var_op == VariationOperator::FirstDeriv) {
             Real const qx = m_lambda * m_gx * m_gx;
             Real const qy = m_lambda * m_gy * m_gy;
             for (int j=m_jlo; j<=m_jhi+1; ++j) {
                 for (int i=m_ilo; i<=m_ihi; ++i) {
                     d_arr(i  ,j,0) += qx;
                     d_arr(i+1,j,0) += qx;
                 }
             }
             for (int j=m_jlo; j<=m_jhi; ++j) {
                 for (int i=m_ilo; i<=m_ihi+1; ++i) {
                     d_arr(i,j  ,0) += qy;
                     d_arr(i,j+1,0) += qy;
                 }
             }
         } else {
             Real const qx = m_lambda * m_gxx * m_gxx;
             Real const qy = m_lambda * m_gyy * m_gyy;
             Real const qc = m_lambda * four * (m_gxx + m_gyy) * (m_gxx + m_gyy);
             for (int j=m_jlo+1; j<=m_jhi; ++j) {
                 for (int i=m_ilo+1; i<=m_ihi; ++i) {
                     d_arr(i-1,j  ,0) += qx;
                     d_arr(i+1,j  ,0) += qx;
                     d_arr(i  ,j-1,0) += qy;
                     d_arr(i  ,j+1,0) += qy;
                     d_arr(i  ,j  ,0) += qc;
                 }
             }
         }
     }

Referenced by conjugateGradient().

Here is the caller graph for this function:

◆ conjugateGradient()

SolveInfo NodalReconstruction::conjugateGradient	(	amrex::FArrayBox const &	rhs,
		amrex::FArrayBox &	x,
		Real	relative_tolerance,
		int	max_iterations
	)		const

inlineprivate

Jacobi-preconditioned conjugate gradient on the SPD operator M. x is both the initial guess and the result; it is reset to zero here because every attempt restarts the correction from the interpolant.

     {
         amrex::FArrayBox diag(m_nodal_box,1,amrex::The_Managed_Arena());
         computeDiagonal(diag);
  
         amrex::FArrayBox r (m_nodal_box,1,amrex::The_Managed_Arena());
         amrex::FArrayBox z (m_nodal_box,1,amrex::The_Managed_Arena());
         amrex::FArrayBox p (m_nodal_box,1,amrex::The_Managed_Arena());
         amrex::FArrayBox Mp(m_nodal_box,1,amrex::The_Managed_Arena());
  
         x.setVal<amrex::RunOn::Host>(zero);
         r.copy<amrex::RunOn::Host>(rhs,0,0,1);   // r = rhs - M*0
         applyJacobi(diag,r,z);
         p.copy<amrex::RunOn::Host>(z,0,0,1);
  
         Real rz = dot(r,z);
         Real const initial = std::sqrt(dot(r,r));
         Real const tol = relative_tolerance * std::max({initial,one});
  
         SolveInfo info;
         info.final_residual = initial;
         if (initial <= tol) {
             info.converged = true;
             return info;
         }
  
         for (int it=0; it<max_iterations; ++it) {
             applyM(p,Mp);
             Real const pMp = dot(p,Mp);
             if (!(pMp > zero) || !std::isfinite(pMp)) {
                 throw std::runtime_error("CG curvature failure.");
             }
  
             Real const alpha = rz/pMp;
             {
                 auto const x_arr  = x.array();
                 auto const r_arr  = r.array();
                 auto const p_arr  = p.const_array();
                 auto const Mp_arr = Mp.const_array();
                 for (int j=m_jlo; j<=m_jhi+1; ++j) {
                     for (int i=m_ilo; i<=m_ihi+1; ++i) {
                         x_arr(i,j,0) += alpha * p_arr(i,j,0);
                         r_arr(i,j,0) -= alpha * Mp_arr(i,j,0);
                     }
                 }
             }
  
             info.iterations = it+1;
             info.final_residual = std::sqrt(dot(r,r));
             if (info.final_residual <= tol) {
                 info.converged = true;
                 return info;
             }
  
             applyJacobi(diag,r,z);
             Real const rz_new = dot(r,z);
             Real const beta = rz_new/rz;
             {
                 auto const p_arr = p.array();
                 auto const z_arr = z.const_array();
                 for (int j=m_jlo; j<=m_jhi+1; ++j) {
                     for (int i=m_ilo; i<=m_ihi+1; ++i) {
                         p_arr(i,j,0) = z_arr(i,j,0) + beta * p_arr(i,j,0);
                     }
                 }
             }
             rz = rz_new;
         }
         return info;
     }

Referenced by solve().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ dot()

Real NodalReconstruction::dot	(	amrex::FArrayBox const &	a,
		amrex::FArrayBox const &	b
	)		const

inlineprivate

     {
         auto const a_arr = a.const_array();
         auto const b_arr = b.const_array();
         Real v = zero;
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 v += a_arr(i,j,0) * b_arr(i,j,0);
             }
         }
         return v;
     }

Referenced by conjugateGradient().

Here is the caller graph for this function:

◆ makeReference()

amrex::FArrayBox NodalReconstruction::makeReference ( amrex::FArrayBox const & T ) const

inline

Direct interpolation of the cell-centered data to the nodes: the average of the (up to four) cells that touch each node, with edge and corner nodes taking the average of the cells that exist. This is the field the solve corrects, and the field it falls back to.

     {
         amrex::FArrayBox R(m_nodal_box,1,amrex::The_Managed_Arena());
         auto const T_arr = T.const_array();
         auto const R_arr = R.array();
  
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             int const jm = std::max(m_jlo, std::min(j-1, m_jhi));
             int const jp = std::max(m_jlo, std::min(j  , m_jhi));
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 int const im = std::max(m_ilo, std::min(i-1, m_ihi));
                 int const ip = std::max(m_ilo, std::min(i  , m_ihi));
                 R_arr(i,j,0) = fourth * ( T_arr(im,jm,0) + T_arr(ip,jm,0)
                                         + T_arr(im,jp,0) + T_arr(ip,jp,0) );
             }
         }
         return R;
     }

◆ maxAbs()

Real NodalReconstruction::maxAbs ( amrex::FArrayBox const & F ) const

inlineprivate

     {
         auto const F_arr = F.const_array();
         Real v = zero;
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 v = std::max(v, std::abs(F_arr(i,j,0)));
             }
         }
         return v;
     }

Referenced by solve().

Here is the caller graph for this function:

◆ maxAverageError()

Real NodalReconstruction::maxAverageError	(	amrex::FArrayBox const &	T,
		amrex::FArrayBox const &	S
	)		const

inline

max_ij | avg4(S) - T |

     {
         auto const T_arr = T.const_array();
         auto const S_arr = S.const_array();
         Real err = zero;
         for (int j=m_jlo; j<=m_jhi; ++j) {
             for (int i=m_ilo; i<=m_ihi; ++i) {
                 Real const avg = fourth * ( S_arr(i  ,j  ,0) + S_arr(i+1,j  ,0)
                                           + S_arr(i  ,j+1,0) + S_arr(i+1,j+1,0) );
                 err = std::max(err, std::abs(avg - T_arr(i,j,0)));
             }
         }
         return err;
     }

Referenced by computeDiagnostics().

Here is the caller graph for this function:

◆ nodalBox()

amrex::Box const& NodalReconstruction::nodalBox ( ) const

inlinenoexcept

                                                              {
         return m_nodal_box;
     }

◆ solve()

std::pair<amrex::FArrayBox,SolveInfo> NodalReconstruction::solve	(	amrex::FArrayBox const &	T,
		amrex::FArrayBox const &	reference,
		VariationOperator	var_op,
		Real	relative_tolerance = `Real(1.e-10)`,
		Real	smoothing = `Real(1.e-2)`,
		Real	min_regularization = `Real(2.5e-4)`,
		int	max_iterations = `1000`
	)

inline

Solve the regularized least-squares problem described above.

Parameters

T	cell-centered data to be matched
reference	the interpolant to correct (see makeReference)
var_op	roughness operator used as the smoothness prior
relative_tolerance	CG stopping tolerance
smoothing	lambda; weight of the roughness penalty
min_regularization	smallest mu tried; raised by factors of four until the correction respects its cap
max_iterations	CG iteration cap per attempt

     {
         if (!(smoothing > zero)) {
             throw std::invalid_argument("Smoothing weight must be positive.");
         }
         if (!(min_regularization > zero)) {
             // With mu == 0 the operator can be singular: the roughness penalty
             // alone does not control every mode (the Laplacian penalty, for
             // instance, is blind to bilinear and discrete-harmonic fields).
             throw std::invalid_argument("Regularization must be positive.");
         }
  
         m_var_op = var_op;
         m_lambda = smoothing;
  
         m_cc_scratch.resize(m_face_box ,1,amrex::The_Managed_Arena());
         m_nd_scratch.resize(m_nodal_box,1,amrex::The_Managed_Arena());
  
         SolveInfo info;
  
         // Misfit of the plain interpolation, and the cap it earns the solve.
         amrex::FArrayBox misfit(m_face_box,1,amrex::The_Managed_Arena());
         Real misfit_max = zero;
         Real t_min      =  std::numeric_limits<Real>::max();
         Real t_max      = -std::numeric_limits<Real>::max();
         {
             applyA(reference,misfit);
             auto const T_arr = T.const_array();
             auto const m_arr = misfit.array();
             for (int j=m_jlo; j<=m_jhi; ++j) {
                 for (int i=m_ilo; i<=m_ihi; ++i) {
                     m_arr(i,j,0) = T_arr(i,j,0) - m_arr(i,j,0);
                     misfit_max = std::max(misfit_max, std::abs(m_arr(i,j,0)));
                     t_min      = std::min(t_min, T_arr(i,j,0));
                     t_max      = std::max(t_max, T_arr(i,j,0));
                 }
             }
         }
         info.interp_average_error = misfit_max;
         info.deviation_cap = std::min(deviation_factor * misfit_max,
                                       relief_factor * (t_max - t_min));
  
         // rhs = A^T (T - A S_ref); the correction solves M E = rhs from E = 0.
         amrex::FArrayBox rhs(m_nodal_box,1,amrex::The_Managed_Arena());
         applyAT(misfit,rhs);
  
         amrex::FArrayBox E(m_nodal_box,1,amrex::The_Managed_Arena());
         amrex::FArrayBox S(m_nodal_box,1,amrex::The_Managed_Arena());
  
         m_mu = min_regularization;
         bool accepted = false;
         for (int attempt=0; attempt<max_attempts; ++attempt) {
             SolveInfo attempt_info = conjugateGradient(rhs,E,relative_tolerance,max_iterations);
             info.iterations += attempt_info.iterations;
             info.final_residual = attempt_info.final_residual;
             info.converged      = attempt_info.converged;
  
             info.deviation = maxAbs(E);
             if (info.deviation <= info.deviation_cap) { accepted = true; break; }
  
             m_mu *= mu_growth;
             ++info.refinements;
         }
  
         if (accepted) {
             auto const S_arr = S.array();
             auto const r_arr = reference.const_array();
             auto const E_arr = E.const_array();
             for (int j=m_jlo; j<=m_jhi+1; ++j) {
                 for (int i=m_ilo; i<=m_ihi+1; ++i) {
                     S_arr(i,j,0) = r_arr(i,j,0) + E_arr(i,j,0);
                 }
             }
         } else {
             // Should not happen -- ||E|| -> 0 as mu -> infinity -- but if it
             // does, the interpolant is a safe, bounded answer.
             S.copy<amrex::RunOn::Host>(reference,0,0,1);
             info.deviation = zero;
         }
         info.regularization = m_mu;
  
         computeDiagnostics(T,S,info);
         return {std::move(S),info};
     }

Here is the call graph for this function:

◆ totalSquaredVariation()

Real NodalReconstruction::totalSquaredVariation ( amrex::FArrayBox const & S ) const

inline

     {
         auto const S_arr = S.const_array();
         Real value = zero;
  
         for (int j=m_jlo; j<=m_jhi+1; ++j) {
             for (int i=m_ilo; i<=m_ihi; ++i) {
                 Real const d = S_arr(i+1,j,0) - S_arr(i,j,0);
                 value += d*d;
             }
         }
         for (int j=m_jlo; j<=m_jhi; ++j) {
             for (int i=m_ilo; i<=m_ihi+1; ++i) {
                 Real const d = S_arr(i,j+1,0) - S_arr(i,j,0);
                 value += d*d;
             }
         }
         return value;
     }

Member Data Documentation

◆ deviation_factor

constexpr amrex::Real NodalReconstruction::deviation_factor = amrex::Real(2.0)

staticconstexprprivate

The correction may not move a node further from the direct interpolation than this multiple of the interpolation misfit ...

Referenced by solve().

◆ m_cc_scratch

amrex::FArrayBox NodalReconstruction::m_cc_scratch

mutableprivate

Referenced by applyM(), and solve().

◆ m_face_box

amrex::Box NodalReconstruction::m_face_box

private

Referenced by solve().

◆ m_gx

amrex::Real NodalReconstruction::m_gx = one

private

Referenced by applyD(), applyDT(), computeDiagonal(), and NodalReconstruction().

◆ m_gxx

amrex::Real NodalReconstruction::m_gxx = one

private

Referenced by applyL(), applyLT(), computeDiagonal(), and NodalReconstruction().

◆ m_gy

amrex::Real NodalReconstruction::m_gy = one

private

Referenced by applyD(), applyDT(), computeDiagonal(), and NodalReconstruction().

◆ m_gyy

amrex::Real NodalReconstruction::m_gyy = one

private

Referenced by applyL(), applyLT(), computeDiagonal(), and NodalReconstruction().

◆ m_ihi

int NodalReconstruction::m_ihi =-1

private

Referenced by applyA(), applyAT(), applyD(), applyDT(), applyJacobi(), applyL(), applyLT(), applyM(), computeDiagnostics(), computeDiagonal(), conjugateGradient(), dot(), makeReference(), maxAbs(), maxAverageError(), solve(), and totalSquaredVariation().

◆ m_ilo

int NodalReconstruction::m_ilo =0

private

Referenced by applyA(), applyAT(), applyD(), applyDT(), applyJacobi(), applyL(), applyLT(), applyM(), computeDiagnostics(), computeDiagonal(), conjugateGradient(), dot(), makeReference(), maxAbs(), maxAverageError(), solve(), and totalSquaredVariation().

◆ m_jhi

int NodalReconstruction::m_jhi =-1

private

Referenced by applyA(), applyAT(), applyD(), applyDT(), applyJacobi(), applyL(), applyLT(), applyM(), computeDiagnostics(), computeDiagonal(), conjugateGradient(), dot(), makeReference(), maxAbs(), maxAverageError(), solve(), and totalSquaredVariation().

◆ m_jlo

int NodalReconstruction::m_jlo =0

private

Referenced by applyA(), applyAT(), applyD(), applyDT(), applyJacobi(), applyL(), applyLT(), applyM(), computeDiagnostics(), computeDiagonal(), conjugateGradient(), dot(), makeReference(), maxAbs(), maxAverageError(), solve(), and totalSquaredVariation().

◆ m_lambda

amrex::Real NodalReconstruction::m_lambda = zero

private

Referenced by applyM(), computeDiagonal(), and solve().

◆ m_mu

amrex::Real NodalReconstruction::m_mu = zero

private

Referenced by applyM(), computeDiagonal(), and solve().

◆ m_nd_scratch

amrex::FArrayBox NodalReconstruction::m_nd_scratch

mutableprivate

Referenced by applyM(), and solve().

◆ m_nodal_box

amrex::Box NodalReconstruction::m_nodal_box

private

Referenced by applyDTD(), applyLTL(), conjugateGradient(), makeReference(), nodalBox(), NodalReconstruction(), and solve().

◆ m_var_op

VariationOperator NodalReconstruction::m_var_op = VariationOperator::FirstDeriv

private

Referenced by applyVariationOperator(), computeDiagonal(), and solve().

◆ max_attempts

constexpr int NodalReconstruction::max_attempts = 24

staticconstexprprivate

Referenced by solve().

◆ mu_growth

constexpr amrex::Real NodalReconstruction::mu_growth = amrex::Real(4.0)

staticconstexprprivate

Factor by which mu is raised when the cap is violated, and the number of times we are willing to do so.

Referenced by solve().

◆ relief_factor

constexpr amrex::Real NodalReconstruction::relief_factor = amrex::Real(0.25)

staticconstexprprivate

... nor than this fraction of the relief of the data.

Referenced by solve().

The documentation for this class was generated from the following file:

Source/Utils/ERF_NodalReconstruction.H

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Why we do not simply invert the averaging operator

What we do instead

Choosing mu

Member Typedef Documentation

◆ Real

Constructor & Destructor Documentation

◆ NodalReconstruction()

Member Function Documentation

◆ applyA()

◆ applyAT()

◆ applyD()

◆ applyDT()

◆ applyDTD()

◆ applyJacobi()

◆ applyL()

◆ applyLT()

◆ applyLTL()

◆ applyM()

◆ applyVariationOperator()

◆ computeDiagnostics()

◆ computeDiagonal()

◆ conjugateGradient()

◆ dot()

◆ makeReference()

◆ maxAbs()

◆ maxAverageError()

◆ nodalBox()

◆ solve()

◆ totalSquaredVariation()

Member Data Documentation

◆ deviation_factor

◆ m_cc_scratch

◆ m_face_box

◆ m_gx

◆ m_gxx

◆ m_gy

◆ m_gyy

◆ m_ihi

◆ m_ilo

◆ m_jhi

◆ m_jlo

◆ m_lambda

◆ m_mu

◆ m_nd_scratch

◆ m_nodal_box

◆ m_var_op

◆ max_attempts

◆ mu_growth

◆ relief_factor