adjoint_curl_curl Module Reference

Adjoint of the curl-curl term in the Pn/Pn pressure residual. More...

Functions/Subroutines
subroutine, public	adjoint_curl_curl_check (prim_sym, adj_sym, cyclic)
	Stop with an error if the boundary conditions make the adjoint curl-curl load differ from the transpose of the primal's curl-curl pressure term.
subroutine, public	adjoint_curl_curl_load (load_x, load_y, load_z, p_adj, coef, gs, mu, rho, work1, work2, work3, work4)
	Compute the adjoint curl-curl load \( \mathbf{L} = \frac{\mu}{\rho} D_c^T M^T D_c^T M^T G p^\dagger \).

Detailed Description

The Karniadakis-Israeli-Orszag splitting of Neko's \( P_N-P_N \) scheme puts the rotational form of the viscous term, evaluated on the extrapolated velocity \( \mathbf{u}_e \), into the pressure residual (pnpn_prs_res). Before gather-scatter, that contribution is

\[ r_{cc}(\mathbf{u}_e) = -\frac{\mu}{\rho} G^T M \mathbf{W}, \qquad \mathbf{W} = M D_c M D_c \mathbf{u}_e, \]

where \( D_c \) is the element-local, pointwise curl, \( M \mathbf{y} = \bar{B}^{-1} \mathrm{gs}(B \mathbf{y}) \) is the mass-weighted average that Neko's curl applies to it (so that curl is \( M D_c \)), with \( B \) the element-local diagonal mass matrix and \( \bar{B} = \mathrm{gs}(B) \) the assembled one (coefBinv holds \( \bar{B}^{-1} \)), and \( G \) is the weak gradient opgrad, \( G_i p = B D_i p \), whose transpose is cdtp. Since \( \bar{B} = \mathrm{gs}(B) \), \( M \) maps every field to a continuous one and leaves continuous fields unchanged, so it is a projection, \( M^2 = M \), and

\[ r_{cc}(\mathbf{u}_e) = -\frac{\mu}{\rho} G^T M D_c M D_c \mathbf{u}_e . \]

The adjoint of this term with respect to \( \mathbf{u}_e \), acting on the adjoint pressure \( p^\dagger \), therefore carries two factors \( M^T \) and is the volume load

\[ \mathbf{L} = \frac{\mu}{\rho} D_c^T M^T D_c^T M^T G p^\dagger, \qquad M^T \mathbf{y} = B \, \mathrm{gs}(\bar{B}^{-1} \mathbf{y}), \]

i.e. for every continuous \( \mathbf{u}_e \) and \( p^\dagger \), \( \mathbf{u}_e \cdot \mathbf{L} = -p^\dagger \cdot r_{cc}(\mathbf{u}_e) \), with both dot products summed over all local (non-unique) nodes and neither \( \mathbf{L} \) nor \( r_{cc} \) gather-scattered.

The transpose is exact in three dimensions only. For a two-dimensional mesh Neko's curl drops the \( z \) derivatives and does not mass-average the \( z \) component, whereas \( \mathbf{L} \) applies the full three-dimensional \( D_c^T \) and \( M^T \) to all three components. The adjoint case therefore rejects two-dimensional meshes, and adjoint_curl_curl_check rejects the boundary conditions for which the load is not the transpose of the primal term.

In the continuous setting \( \nabla \times \nabla p^\dagger = 0 \) and the term reduces to a boundary integral. The discrete operators do not satisfy that identity, so a boundary-only form is not the transpose of what the primal computes; the load acts on the whole volume.

Function/Subroutine Documentation

adjoint_curl_curl_check()

subroutine, public adjoint_curl_curl::adjoint_curl_curl_check	(	type(facet_normal_t), intent(in)	prim_sym,
		type(facet_normal_t), intent(in)	adj_sym,
		logical, intent(in)	cyclic
	)

The load does not transpose the symmetry-surface part of the primal term (bc_sym_surface), nor the cyclic rotations applied inside the primal's curl. The symmetry test is reduced over all ranks, so that every rank stops, also those that hold no symmetry facet. The restriction to three-dimensional meshes (see the module description) is not tested here: it has to be enforced before the adjoint fluid is initialised, which assumes a three-dimensional mesh.

Parameters

prim_sym	The symmetry surface of the primal pressure residual.
adj_sym	The symmetry surface of the adjoint pressure residual.
cyclic	Whether the primal has cyclic periodic boundaries.

Definition at line 118 of file adjoint_curl_curl.f90.

adjoint_curl_curl_load()

subroutine, public adjoint_curl_curl::adjoint_curl_curl_load	(	type(field_t), intent(inout)	load_x,
		type(field_t), intent(inout)	load_y,
		type(field_t), intent(inout)	load_z,
		type(field_t), intent(in)	p_adj,
		type(coef_t), intent(in)	coef,
		type(gs_t), intent(inout)	gs,
		real(kind=rp), intent(in)	mu,
		real(kind=rp), intent(in)	rho,
		type(field_t), intent(inout)	work1,
		type(field_t), intent(inout)	work2,
		type(field_t), intent(inout)	work3,
		type(field_t), intent(inout)	work4
	)

See the module description for the notation. The load is element-local (not gather-scattered), since it is added to the adjoint forcing, which the velocity residual gather-scatters. It is added to the forcing before makeabf, which scales the forcing by \( \rho \) and extrapolates it in time, so the velocity right-hand side receives \( \mu D_c^T M^T D_c^T M^T G p^\dagger \), weighted over the lagged adjoint pressures as the primal's extrapolation weights \( \mathbf{u}_e \). The material properties are taken as constant, as in the primal residual. The cyclic rotations and the symmetry-surface part of the primal term are not transposed, so the load is not valid for cyclic or symmetry boundaries (see adjoint_curl_curl_check), nor for a two-dimensional mesh. All field arguments must be distinct.

Parameters

load_x	x component of the load.
load_y	y component of the load.
load_z	z component of the load.
p_adj	The adjoint pressure (continuous).
coef	The SEM coefficients.
gs	The gather-scatter operator of the velocity space.
mu	The constant dynamic viscosity.
rho	The constant density.
work1	Work field.
work2	Work field.
work3	Work field.
work4	Work field.

Definition at line 170 of file adjoint_curl_curl.f90.