SELF_ColumnOperators_3D_t Module

  1. SELF_ColumnOperators_3D_t.f90
  2. SELF_ColumnOperators_3D_t

Column operators on extruded meshes (Mesh3D%ExtrudedMesh).

An extruded mesh stacks nz = mesh%nLayers hexahedra above every surface element e2 of a Mesh2D. Local element e3 = (e2-1)*nz + k is layer k of that column, k = 1 at the bottom, and a column is never split across MPI ranks. Every element maps

x(xi,eta,zeta) = ( X(xi,eta), Y(xi,eta), Z(xi,eta,zeta) ),

so the lines of constant (xi,eta), the third index of SELF's interior arrays, are vertical. The operators below act along those lines. They are local to a column, so none of them communicates.

Notation. N is the polynomial degree, w_l (l = 0..N) the quadrature weights, M = diag(w), D the nodal derivative matrix, D(row=i,col=m) = interp%dMatrix(m,i), e_0 and e_N the unit vectors of the bottom (zeta = -1) and the top (zeta = 1) node, and Z_zeta = dz/dzeta, the (z,zeta) entry dxds%interior(...,1,3,3) of the covariant basis SEMHex stores at the solution nodes. A superscript (k) is the layer.

ColumnIntegral, volume to surface, in the unit of q times metres:

(C q)ij = sum w_l (Z_zeta)^(k)_ijl q^(k)_ijl .} sum_{l=0..N

Broadcast, surface to volume, unit unchanged:

(B Q)^(k)_ijl = Q_ij .

Because the Jacobian of an extruded element factors as J = J_2 Z_zeta at every node, with J_2 the surface Jacobian, B is the adjoint of C,

_3 = _2 , _3 = sum J w_i w_j w_l a c , _2 = sum J_2 w_i w_j A C .

VerticalSweepUp and VerticalSweepDown solve, on every vertical line of every element, the one dimensional upwind DG problems of SELF-HPE's continuity equation and hydrostatic relation. Upward, from a given bottom trace,

A q^(k) = r^(k) + M^{-1} e_0 q_in^(k) , A = D + M^{-1} e_0 e_0^T , q_in^(1) = qBottom , q_in^(k) = q^(k-1)_N ;

downward, from a given top trace,

At q^(k) = r^(k) - M^{-1} e_N q_in^(k) , At = D - M^{-1} e_N e_N^T , q_in^(nz) = qTop , q_in^(k) = q^(k+1)_0 .

The right hand side r is in reference scaling: the sweeps solve dq/dzeta = r in the reference coordinate of each element, and they apply no metric factor. A caller that solves a physical equation such as J2c dw/dzeta = -(horizontal terms) folds the metric factors into r before the call. Since A 1 = M^{-1} e_0 and At 1 = -M^{-1} e_N, the solution is the inflow value plus an element-local particular solution,

q^(k) = q_in^(k) 1 + A^{-1} r^(k) (upward), q^(k) = q_in^(k) 1 + At^{-1} r^(k) (downward),

with A^{-1} and At^{-1} precomputed once in Init. The layers of a column are visited in order, because each takes its inflow from the one before; columns, the nodes (i,j) of a column and the variables are independent. By the summation-by-parts property M D + D^T M = e_N e_N^T - e_0 e_0^T of Gauss-Lobatto nodes the two sweeps are adjoint: with w from the upward sweep and p from the downward sweep, on every vertical line,

sum_k [ p^(k)T M r_w^(k) + w^(k)T M r_p^(k) ] = p_in^(nz) w^(nz)_N - p^(1)_0 w_in^(1) .

The sweeps need nodes at zeta = -1 and zeta = 1, so they require Gauss-Lobatto control nodes. The column integral and the broadcast work with any quadrature.

Fields. Volume fields are Scalar3D, or any extension such as MappedScalar3D, on the extruded mesh with nElem = nz*nColumns. Surface fields are Scalar2D, or any extension such as MappedScalar2D, with nElem = nColumns, the same interpolant degree and the same number of variables; surface element e2 is the column above the caller's surface element e2. Only nodal values are touched, so neither field needs a geometry. On GPU builds every operator works on the device copies (interior_gpu) and leaves the host copies untouched.


Uses

  • module~~self_columnoperators_3d_t~~UsesGraph module~self_columnoperators_3d_t SELF_ColumnOperators_3D_t module~self_geometry_3d SELF_Geometry_3D module~self_columnoperators_3d_t->module~self_geometry_3d module~self_scalar_2d SELF_Scalar_2D module~self_columnoperators_3d_t->module~self_scalar_2d module~self_scalar_3d SELF_Scalar_3D module~self_columnoperators_3d_t->module~self_scalar_3d iso_fortran_env iso_fortran_env module~self_columnoperators_3d_t->iso_fortran_env module~self_constants SELF_Constants module~self_columnoperators_3d_t->module~self_constants module~self_mesh_3d SELF_Mesh_3D module~self_columnoperators_3d_t->module~self_mesh_3d module~self_lagrange~3 SELF_Lagrange module~self_columnoperators_3d_t->module~self_lagrange~3 module~self_geometry_3d->module~self_scalar_3d module~self_geometry_3d->module~self_constants module~self_geometry_3d->module~self_mesh_3d module~self_geometry_3d->module~self_lagrange~3 module~self_vector_3d SELF_Vector_3D module~self_geometry_3d->module~self_vector_3d module~self_data SELF_Data module~self_geometry_3d->module~self_data module~self_supportroutines SELF_SupportRoutines module~self_geometry_3d->module~self_supportroutines module~self_tensor_3d SELF_Tensor_3D module~self_geometry_3d->module~self_tensor_3d module~self_scalar_2d->module~self_constants module~self_gpuinterfaces SELF_GPUInterfaces module~self_scalar_2d->module~self_gpuinterfaces module~self_scalar_2d_t SELF_Scalar_2D_t module~self_scalar_2d->module~self_scalar_2d_t iso_c_binding iso_c_binding module~self_scalar_2d->iso_c_binding module~self_gpu~2 SELF_GPU module~self_scalar_2d->module~self_gpu~2 module~self_scalar_3d->module~self_constants module~self_scalar_3d->module~self_gpuinterfaces module~self_scalar_3d->iso_c_binding module~self_scalar_3d->module~self_gpu~2 module~self_scalar_3d_t SELF_Scalar_3D_t module~self_scalar_3d->module~self_scalar_3d_t module~self_constants->iso_fortran_env module~self_constants->iso_c_binding module~self_mesh_3d->iso_c_binding module~self_mesh_3d->module~self_gpu~2 module~self_mesh_3d_t SELF_Mesh_3D_t module~self_mesh_3d->module~self_mesh_3d_t module~self_lagrange~3->iso_fortran_env module~self_lagrange~3->module~self_constants module~self_lagrange_t SELF_Lagrange_t module~self_lagrange~3->module~self_lagrange_t module~self_lagrange~3->iso_c_binding module~self_lagrange~3->module~self_gpu~2 module~self_lagrange_t->iso_fortran_env module~self_lagrange_t->module~self_constants module~self_lagrange_t->module~self_supportroutines module~self_lagrange_t->iso_c_binding HDF5 HDF5 module~self_lagrange_t->HDF5 module~self_hdf5 SELF_HDF5 module~self_lagrange_t->module~self_hdf5 module~self_quadrature SELF_Quadrature module~self_lagrange_t->module~self_quadrature module~self_vector_3d->module~self_constants module~self_vector_3d->module~self_gpuinterfaces module~self_vector_3d->iso_c_binding module~self_vector_3d->module~self_gpu~2 module~self_vector_3d_t SELF_Vector_3D_t module~self_vector_3d->module~self_vector_3d_t module~self_data->module~self_constants module~self_data->module~self_lagrange~3 module~self_data->iso_c_binding module~self_data->HDF5 module~self_metadata SELF_Metadata module~self_data->module~self_metadata module~self_data->module~self_hdf5 FEQParse FEQParse module~self_data->FEQParse module~self_supportroutines->iso_fortran_env module~self_supportroutines->module~self_constants module~self_gpuinterfaces->iso_c_binding module~self_gpuinterfaces->module~self_gpu~2 module~self_tensor_3d->iso_fortran_env module~self_tensor_3d->module~self_constants module~self_tensor_3d->iso_c_binding module~self_tensor_3d->module~self_gpu~2 module~self_tensor_3d_t SELF_Tensor_3D_t module~self_tensor_3d->module~self_tensor_3d_t module~self_scalar_2d_t->module~self_constants module~self_scalar_2d_t->module~self_lagrange~3 module~self_scalar_2d_t->module~self_data module~self_scalar_2d_t->iso_c_binding module~self_scalar_2d_t->HDF5 module~self_scalar_2d_t->module~self_metadata module~self_scalar_2d_t->module~self_hdf5 module~self_scalar_2d_t->FEQParse module~self_datapool SELF_DataPool module~self_scalar_2d_t->module~self_datapool module~self_gpu~2->iso_c_binding module~self_gpu_enums~2 SELF_GPU_enums module~self_gpu~2->module~self_gpu_enums~2 module~self_scalar_3d_t->module~self_constants module~self_scalar_3d_t->module~self_lagrange~3 module~self_scalar_3d_t->module~self_data module~self_scalar_3d_t->iso_c_binding module~self_scalar_3d_t->HDF5 module~self_scalar_3d_t->module~self_metadata module~self_scalar_3d_t->module~self_hdf5 module~self_scalar_3d_t->FEQParse module~self_scalar_3d_t->module~self_datapool module~self_mesh_3d_t->module~self_constants module~self_mesh_3d_t->module~self_lagrange~3 module~self_mesh_3d_t->module~self_supportroutines module~self_mesh_3d_t->iso_c_binding module~self_mesh_3d_t->HDF5 module~self_mesh_3d_t->module~self_hdf5 module~self_mesh_3d_t->module~self_quadrature module~self_mesh_2d_t SELF_Mesh_2D_t module~self_mesh_3d_t->module~self_mesh_2d_t module~self_domaindecomposition~2 SELF_DomainDecomposition module~self_mesh_3d_t->module~self_domaindecomposition~2 module~self_mesh SELF_Mesh module~self_mesh_3d_t->module~self_mesh module~self_metadata->HDF5 module~self_metadata->module~self_hdf5 module~self_hdf5->iso_fortran_env module~self_hdf5->module~self_constants module~self_hdf5->HDF5 mpi mpi module~self_hdf5->mpi module~self_vector_3d_t->module~self_constants module~self_vector_3d_t->module~self_lagrange~3 module~self_vector_3d_t->module~self_data module~self_vector_3d_t->iso_c_binding module~self_vector_3d_t->HDF5 module~self_vector_3d_t->module~self_metadata module~self_vector_3d_t->module~self_hdf5 module~self_vector_3d_t->FEQParse module~self_vector_3d_t->module~self_datapool module~self_quadrature->iso_fortran_env module~self_quadrature->module~self_constants module~self_tensor_3d_t->module~self_constants module~self_tensor_3d_t->module~self_lagrange~3 module~self_tensor_3d_t->module~self_data module~self_tensor_3d_t->iso_c_binding module~self_tensor_3d_t->HDF5 module~self_tensor_3d_t->module~self_metadata module~self_tensor_3d_t->module~self_hdf5 module~self_tensor_3d_t->FEQParse module~self_tensor_3d_t->module~self_datapool module~self_datapool->module~self_constants module~self_gpu_enums~2->iso_c_binding module~self_mesh_2d_t->module~self_constants module~self_mesh_2d_t->module~self_lagrange~3 module~self_mesh_2d_t->module~self_supportroutines module~self_mesh_2d_t->iso_c_binding module~self_mesh_2d_t->HDF5 module~self_mesh_2d_t->module~self_hdf5 module~self_mesh_2d_t->module~self_quadrature module~self_mesh_2d_t->module~self_domaindecomposition~2 module~self_mesh_2d_t->module~self_mesh module~self_domaindecomposition~2->iso_fortran_env module~self_domaindecomposition~2->iso_c_binding module~self_domaindecomposition~2->mpi module~self_domaindecomposition_t SELF_DomainDecomposition_t module~self_domaindecomposition~2->module~self_domaindecomposition_t module~self_mesh->module~self_constants module~self_mesh->iso_c_binding module~self_mesh->module~self_domaindecomposition~2 module~self_domaindecomposition_t->module~self_constants module~self_domaindecomposition_t->module~self_lagrange~3 module~self_domaindecomposition_t->module~self_supportroutines module~self_domaindecomposition_t->iso_c_binding module~self_domaindecomposition_t->mpi

Used by

  • module~~self_columnoperators_3d_t~~UsedByGraph module~self_columnoperators_3d_t SELF_ColumnOperators_3D_t module~self_columnoperators_3d SELF_ColumnOperators_3D module~self_columnoperators_3d->module~self_columnoperators_3d_t module~self_columnoperators_3d~2 SELF_ColumnOperators_3D module~self_columnoperators_3d~2->module~self_columnoperators_3d_t

Contents

Derived Types

Subroutines


Derived Types

type, public :: ColumnOperators3D_t

Components

TypeVisibilityAttributesNameInitial
integer, public :: N =0

Polynomial degree of the solution
real(kind=prec), public, pointer, contiguous, dimension(:,:):: downInverse=> null()

At^{-1} in the same convention
type(SEMHex), public, pointer:: geometry=> null()
type(Lagrange), public, pointer:: interp=> null()
logical, public :: lobatto =.false.

True when the control nodes are Gauss-Lobatto, which the sweeps require
integer, public :: nColumns =0

Rank-local columns, which is the rank-local element count of the surface mesh
integer, public :: nElem =0

Rank-local 3-D elements, nLayers*nColumns
integer, public :: nLayers =0

Layers per column, nz
real(kind=prec), public, pointer, contiguous, dimension(:,:):: upInverse=> null()

A^{-1} in SELF's matrix convention: (A^{-1} r)_l = sum_m upInverse(m,l) r_m

Type-Bound Procedures

procedure, public :: Broadcast => Broadcast_ColumnOperators3D_t
procedure, public :: ColumnIntegral => ColumnIntegral_ColumnOperators3D_t
procedure, public :: Free => Free_ColumnOperators3D_t
procedure, public :: Init => Init_ColumnOperators3D_t
procedure, public :: VerticalSweepDown => VerticalSweepDown_ColumnOperators3D_t
procedure, public :: VerticalSweepUp => VerticalSweepUp_ColumnOperators3D_t

Subroutines

public subroutine Broadcast_ColumnOperators3D_t(this, qbar, q)

Broadcast of a surface field to every node of its column,

Read more…

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(in) :: this
class(Scalar2D), intent(in) :: qbar
class(Scalar3D), intent(inout) :: q

public subroutine CheckFields_ColumnOperators3D_t(this, volume, surface, caller)

Stops when a volume and a surface field do not match the extruded mesh the operators were initialized for, or each other.

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(in) :: this
class(Scalar3D), intent(in) :: volume
class(Scalar2D), intent(in) :: surface
character, intent(in) :: caller

public subroutine CheckSweep_ColumnOperators3D_t(this, r, inflow, q, caller)

The checks of CheckFields for a vertical sweep, plus the requirement of Gauss-Lobatto nodes and of a right hand side and solution of the same shape.

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(in) :: this
class(Scalar3D), intent(in) :: r
class(Scalar2D), intent(in) :: inflow
class(Scalar3D), intent(in) :: q
character, intent(in) :: caller

public subroutine ColumnIntegral_ColumnOperators3D_t(this, q, qbar)

Column integral of a volume field,

Read more…

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(in) :: this
class(Scalar3D), intent(in) :: q
class(Scalar2D), intent(inout) :: qbar

public subroutine Free_ColumnOperators3D_t(this)

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(inout) :: this

public subroutine Init_ColumnOperators3D_t(this, interp, mesh, geometry)

Binds the column operators to an extruded mesh and its geometry and, on Gauss-Lobatto nodes, precomputes the inverses of the two sweep matrices

Read more…

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(out) :: this
type(Lagrange), intent(in), target:: interp
type(Mesh3D), intent(in) :: mesh
type(SEMHex), intent(in), target:: geometry

public subroutine InvertMatrix(a, aInv, n)

Inverts the n x n matrix a by Gauss-Jordan elimination with partial pivoting, in double precision. Used once in Init for the (N+1) x (N+1) sweep matrices, which are small and well enough conditioned for a direct solve; this follows the routine of the same name in SELF_RefinementIndicator_3D_t.

Arguments

TypeIntentOptionalAttributesName
real(kind=real64), intent(in) :: a(1:n,1:n)
real(kind=real64), intent(out) :: aInv(1:n,1:n)
integer, intent(in) :: n

public subroutine VerticalSweepDown_ColumnOperators3D_t(this, r, qTop, q)

Downward sweep: solves, layer by layer from the top of every column,

Read more…

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(in) :: this
class(Scalar3D), intent(in) :: r
class(Scalar2D), intent(in) :: qTop
class(Scalar3D), intent(inout) :: q

public subroutine VerticalSweepUp_ColumnOperators3D_t(this, r, qBottom, q)

Upward sweep: solves, layer by layer from the bottom of every column,

Read more…

Arguments

TypeIntentOptionalAttributesName
class(ColumnOperators3D_t), intent(in) :: this
class(Scalar3D), intent(in) :: r
class(Scalar2D), intent(in) :: qBottom
class(Scalar3D), intent(inout) :: q