recon1d_plm_wls module reference

Piecewise Linear Method using Weighted Conservative Least Squares 1D reconstruction.

Data Types

plm_wls

PLM reconstruction using Weighted Least Squares constrained to conserve for central cell.

Functions/Subroutines

`init()`	Initialize a 1D PLM reconstruction for n cells.
`reconstruct()`	Calculate a 1D PLM_WLS reconstruction based on h(:) and u(:)
`f()`	Value of PLM_WLS reconstruction at a point in cell k [A].
`dfdx()`	Derivative of PLM_WLS reconstruction at a point in cell k [A].
`average()`	Average between xa and xb for cell k of a 1D PLM reconstruction [A].
`destroy()`	Deallocate the PLM reconstruction.
`check_reconstruction()`	Checks the PLM_WLS reconstruction for consistency.
`ls_error()`	Returns local least squares error for a particular cell.
`unit_tests()`	Runs PLM_WLS reconstruction unit tests and returns True for any fails, False otherwise.

Detailed Description

This implementation of PLM fits the slope using least squares, but retains conservation for the central cell by passing through the central value. Cell-wise reconstructions are NOT limited by neighbours. Thus, this reconstruction does not yield monotonic profiles needed for the general remapping problem.

The algorithm solves the least squares problem of fitting a straight line through the neighboring data. The line is constained to pass through the center cell, \((x_{k}, y_{k})\), so that the construction is conservative. The more general function \(f(x) = a_{k} + b_{k} x\) would not conserve for arbitrary data.

The unknown parameter \(b_{k}\) in the line

\[f(x) = y_{k} + b_{k} ( x - x_{k} ) / h_{k}\]

is fit to neighbors \(x_{k-1}, y_{k-1}\) and \(x_{k+1}, y_{k+1}\).

Denoting \(y'_{k+j} = y_{k+j} - y_{k}\) and \(x'_{k+j} = x_{k+j} - x_{k}\) the local error is

\[\begin{split}\begin{align} e_{k+j} &= b_k \frac{ x_{k+j} - x_{k} }{ h_{k} } + y_{k} - y_{k+j} \\\\ &= b_k \frac{ x'_{k+j} }{ h_{k} } - y'_{k+j} \;\; . \end{align}\end{split}\]

We use volume weighting in the loss

\[G(b) = h_{k-1} e_{k-1}^2 + h_{k+1} e_{k+1}^2 \;\; .\]

When solving for \(b_k\), we solve \(dG/db = 0\) where

\[\begin{split}\begin{align} dG/db &= 2 h_{k-1} e_{k-1} \frac{ de_{k-1} }{db} + 2 h_{k+1} e_{k+1} \frac{ de_{k+1} }{db} \\\\ &= 2 h_{k-1} ( b_k \frac{ x'_{k-1} }{ h_{k} } - \frac{ y'_{k-1} ) x'_{k-1} }{ h_{k} } + 2 h_{k+1} ( b_k \frac{ x'_{k+1} }{ h_{k} } - \frac{ y'_{k+1} ) x'_{k+1} }{ h_{k} } \\\\ &= 4 b_k \frac{ < h x'^2 > }{ h_{k}^2 } - 4 \frac{ < h x' y' > }{ h_{k} } \end{align}\end{split}\]

and where \(< a > = \frac{1}{2} ( a_{k-1} + a_{k+1} )\). Thus

\[b_k = \frac{ h_{k} < h x' y' > }{ < h x'^2 > } \;\; .\]

When evaluating the loss, \(G\), some rearrangement is necessary to reduce truncation errors. Since

\[\begin{split}\begin{align} e_{k+j}^2 &= \left( b \frac{ x'_{k+j} }{ h_{k} } - y'_{k+j} \right)^2 \\\\ &= b^2 \frac{ {x'}_{k+j}^2 }{ h_{k}^2 } - 2 b \frac{ x'_{k+j} y'_{k+j} }{ h_{k} } + {y'}_{k+j}^2 \end{align}\end{split}\]

then

\[\begin{split}\begin{align} G(b) &= 2 < h e^2 > \\\\ &= 2 b^2 \frac{ < h {x'}^2 > }{ h_{k}^2 } - 4 b \frac{ < h x' y' > }{ h_{k} } + 2 < h' {y'}^2 > \;\; . \end{align}\end{split}\]

If we denote the value of b that yields the minimum value as \(b^*\) then

\[G(b^*) = < h {y'}^2 > - \frac{ < h x' y' >^2 }{ < h {x'}^2 > } \;\; .\]

Let

\[\begin{split}\begin{align} G''(b) &= G(b) - G(b^*) \\\\ &= b^2 \frac{ < h {x'}^2 > }{ h_{k}^2 } - 2 b \frac{ < h x' y' > }{ h_{k} } + \frac{ < h x' y' > }{ < h {x'}^2 > } \\\\ &= \frac{ \left( b < h {x'}^2 > - h_{k} < h x' y' > \right)^2 }{ h_{k} < h {x'}^2 > } \;\; . \end{align}\end{split}\]

Minimizing \(G''(b)\) is equivalent to minimizing \(G(b)\) for the same data. \(G''(b^*)=0\) so evaluation with the last form, in the vicinity of \(b^*\), avoids large cancelling terms.

Type Documentation

type recon1d_plm_wls/plm_wls

PLM reconstruction using Weighted Least Squares constrained to conserve for central cell.

Type fields:

% ul :: real, dimension(:), allocatable, private Left edge value [A].
% ur :: real, dimension(:), allocatable, private Right edge value [A].
% slp :: real, dimension(:), allocatable, private, private Difference across cell, ur - ul [A]. This is redundant with ul and ur and not used in any evaluations, but is needed for testing.
% init :: procedure, private Implementation of the PLM_WLS initialization.
% reconstruct :: procedure, private Implementation of the PLM_WLS reconstruction.
% average :: procedure, private Implementation of the PLM_WLS average over an interval [A].
% f :: procedure, private Implementation of evaluating the PLM_WLS reconstruction at a point [A].
% dfdx :: procedure, private Implementation of the derivative of the PLM_WLS reconstruction at a point [A].
% destroy :: procedure, private Implementation of deallocation for PLM_WLS.
% check_reconstruction :: procedure, private Implementation of check reconstruction for the PLM_WLS reconstruction.
% unit_tests :: procedure, private Implementation of unit tests for the PLM_WLS reconstruction.
% init_parent :: procedure, private Duplicate interface to
% reconstruct_parent :: procedure, private Duplicate interface to

[source]

Function/Subroutine Documentation

subroutine recon1d_plm_wls/init(this, n, h_neglect, check)

Initialize a 1D PLM reconstruction for n cells.

Parameters:

this :: this [out] This reconstruction
n :: n [in] Number of cells in this column
h_neglect :: h_neglect [in] A negligibly small width used in cell reconstructions [H]
check :: check [in] If true, enable some consistency checking

[source]

subroutine recon1d_plm_wls/reconstruct(this, h, u)

Calculate a 1D PLM_WLS reconstruction based on h(:) and u(:)

Parameters:

this :: this [inout] This reconstruction
h :: h [in] Grid spacing (thickness) [typically H]
u :: u [in] Cell mean values [A]

[source]

function recon1d_plm_wls/f(this, k, x)

Value of PLM_WLS reconstruction at a point in cell k [A].

Parameters:

this :: this [in] This reconstruction
k :: k [in] Cell number
x :: x [in] Non-dimensional position within element [nondim]

[source]

function recon1d_plm_wls/dfdx(this, k, x)

Derivative of PLM_WLS reconstruction at a point in cell k [A].

Parameters:

this :: this [in] This reconstruction
k :: k [in] Cell number
x :: x [in] Non-dimensional position within element [nondim]

[source]

function recon1d_plm_wls/average(this, k, xa, xb)

Average between xa and xb for cell k of a 1D PLM reconstruction [A].

Parameters:

this :: this [in] This reconstruction
k :: k [in] Cell number
xa :: xa [in] Start of averaging interval on element (0 to 1)
xb :: xb [in] End of averaging interval on element (0 to 1)

[source]

subroutine recon1d_plm_wls/destroy(this)

Deallocate the PLM reconstruction.

Parameters:: this :: this [inout] This reconstruction

[source]

function recon1d_plm_wls/check_reconstruction(this, h, u)

Checks the PLM_WLS reconstruction for consistency.

Parameters:

this :: this [in] This reconstruction
h :: h [in] Grid spacing (thickness) [typically H]
u :: u [in] Cell mean values [A]

Call to:

ls_error

[source]

function recon1d_plm_wls/ls_error(this, k, h, u)

Returns local least squares error for a particular cell.

Note that this is the error relative to the minimum of the loss function so that at the true solution this function returns zero. See module documentation.

Parameters:

this :: this [in] This reconstruction
k :: k [in] Cell number
h :: h [in] Grid spacing (thickness) [typically H]
u :: u [in] Cell mean values [A]

Called from:

check_reconstruction unit_tests

[source]

function recon1d_plm_wls/unit_tests(this, verbose, stdout, stderr)

Runs PLM_WLS reconstruction unit tests and returns True for any fails, False otherwise.

Parameters:

this :: this [inout] This reconstruction
verbose :: verbose [in] True, if verbose
stdout :: stdout [in] I/O channel for stdout
stderr :: stderr [in] I/O channel for stderr

Call to:

ls_error

[source]