Matrix-free / matrix-based epsilon-div-div operator on wedge elements in a spherical shell. More...

#include <epsilon_divdiv_kerngen.hpp>

Public Types
using	SrcVectorType = linalg::VectorQ1Vec< ScalarT, VecDim >

using	DstVectorType = linalg::VectorQ1Vec< ScalarT, VecDim >

using	ScalarType = ScalarT

using	Grid4DDataLocalMatrices = terra::grid::Grid4DDataMatrices< ScalarType, LocalMatrixDim, LocalMatrixDim, 2 >

using	LocalMatrixStorage = linalg::solvers::LocalMatrixStorage< ScalarType, LocalMatrixDim >

using	Team = Kokkos::TeamPolicy<>::member_type

Public Member Functions
	EpsilonDivDivKerngen (const grid::shell::DistributedDomain &domain, const grid::Grid3DDataVec< ScalarT, 3 > &grid, const grid::Grid2DDataScalar< ScalarT > &radii, const grid::Grid4DDataScalar< grid::shell::ShellBoundaryFlag > &mask, const grid::Grid4DDataScalar< ScalarT > &k, BoundaryConditions bcs, bool diagonal, linalg::OperatorApplyMode operator_apply_mode=linalg::OperatorApplyMode::Replace, linalg::OperatorCommunicationMode operator_communication_mode=linalg::OperatorCommunicationMode::CommunicateAdditively, linalg::OperatorStoredMatrixMode operator_stored_matrix_mode=linalg::OperatorStoredMatrixMode::Off)

void	set_operator_apply_and_communication_modes (const linalg::OperatorApplyMode operator_apply_mode, const linalg::OperatorCommunicationMode operator_communication_mode)

void	set_diagonal (bool v)

void	set_boundary_conditions (BoundaryConditions bcs)

const grid::Grid4DDataScalar< ScalarType > &	k_grid_data ()

const grid::shell::DistributedDomain &	get_domain () const

grid::Grid2DDataScalar< ScalarT >	get_radii () const

grid::Grid3DDataVec< ScalarT, 3 >	get_grid ()

bool	has_flag (const int local_subdomain_id, const int x_cell, const int y_cell, const int r_cell, grid::shell::ShellBoundaryFlag flag) const

void	set_stored_matrix_mode (linalg::OperatorStoredMatrixMode operator_stored_matrix_mode, int level_range, grid::Grid4DDataScalar< ScalarType > GCAElements)

linalg::OperatorStoredMatrixMode	get_stored_matrix_mode ()

void	set_local_matrix (const int local_subdomain_id, const int x_cell, const int y_cell, const int r_cell, const int wedge, const dense::Mat< ScalarT, LocalMatrixDim, LocalMatrixDim > &mat) const

dense::Mat< ScalarT, LocalMatrixDim, LocalMatrixDim >	get_local_matrix (const int local_subdomain_id, const int x_cell, const int y_cell, const int r_cell, const int wedge) const

void	apply_impl (const SrcVectorType &src, DstVectorType &dst)
	Apply the operator: dst <- Op(src) (or additive, depending on mode).

void	column_grad_to_sym (const int dim, const double g0, const double g1, const double g2, double &E00, double &E11, double &E22, double &sym01, double &sym02, double &sym12, double &gdd) const
	Convert one gradient column into symmetric-gradient and div contributions.

size_t	team_shmem_size (const int) const
	Team scratch memory size for fast paths.

void	operator() (const Team &team) const
	Legacy generic team operator.

void	assemble_trial_test_vecs (const int wedge, const dense::Vec< ScalarType, VecDim > &quad_point, const ScalarType quad_weight, const ScalarT r_1, const ScalarT r_2, dense::Vec< ScalarT, 3 >(wedge_phy_surf)[3], const dense::Vec< ScalarT, 6 > k_local_hex, const int dimi, const int dimj, dense::Mat< ScalarType, VecDim, VecDim > sym_grad_i, dense::Mat< ScalarType, VecDim, VecDim > sym_grad_j, ScalarType &jdet_keval_quadweight) const

dense::Mat< ScalarT, LocalMatrixDim, LocalMatrixDim >	assemble_local_matrix (const int local_subdomain_id, const int x_cell, const int y_cell, const int r_cell, const int wedge) const
	Assemble one wedge-local 18x18 matrix (slow path / on-demand assembly).

Static Public Attributes
static constexpr int	LocalMatrixDim = 18

Detailed Description

template<typename ScalarT, int VecDim = 3>
class terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >

Matrix-free / matrix-based epsilon-div-div operator on wedge elements in a spherical shell.

This class supports three execution paths:

1) Slow path (local matrix path)

Used if local matrices are stored (full or selective)
Reuses assembled/stored 18x18 local matrices

2) Fast path: Dirichlet/Neumann

Matrix-free
Handles Dirichlet by skipping constrained face-node couplings (and diagonal treatment if requested)
Neumann is naturally included

3) Fast path: Free-slip

Matrix-free
Applies trial and test-side tangential projection on free-slip boundaries
Preserves behavior consistent with the slow path for iterative solves

IMPORTANT DESIGN CHANGE:

The path decision is made on the HOST and cached in kernel_path_. apply_impl() dispatches to a different kernel launch per path. This avoids a runtime branch inside the hot device kernel.

Member Typedef Documentation

◆ DstVectorType

template<typename ScalarT , int VecDim = 3>

using terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::DstVectorType = linalg::VectorQ1Vec< ScalarT, VecDim >

◆ Grid4DDataLocalMatrices

template<typename ScalarT , int VecDim = 3>

using terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::Grid4DDataLocalMatrices = terra::grid::Grid4DDataMatrices< ScalarType, LocalMatrixDim, LocalMatrixDim, 2 >

◆ LocalMatrixStorage

template<typename ScalarT , int VecDim = 3>

using terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::LocalMatrixStorage = linalg::solvers::LocalMatrixStorage< ScalarType, LocalMatrixDim >

◆ ScalarType

template<typename ScalarT , int VecDim = 3>

using terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::ScalarType = ScalarT

◆ SrcVectorType

template<typename ScalarT , int VecDim = 3>

using terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::SrcVectorType = linalg::VectorQ1Vec< ScalarT, VecDim >

◆ Team

template<typename ScalarT , int VecDim = 3>

using terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::Team = Kokkos::TeamPolicy<>::member_type

Constructor & Destructor Documentation

◆ EpsilonDivDivKerngen()

template<typename ScalarT , int VecDim = 3>

terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::EpsilonDivDivKerngen	(	const grid::shell::DistributedDomain &	domain,
		const grid::Grid3DDataVec< ScalarT, 3 > &	grid,
		const grid::Grid2DDataScalar< ScalarT > &	radii,
		const grid::Grid4DDataScalar< grid::shell::ShellBoundaryFlag > &	mask,
		const grid::Grid4DDataScalar< ScalarT > &	k,
		BoundaryConditions	bcs,
		bool	diagonal,
		linalg::OperatorApplyMode	operator_apply_mode = `linalg::OperatorApplyMode::Replace`,
		linalg::OperatorCommunicationMode	operator_communication_mode = `linalg::OperatorCommunicationMode::CommunicateAdditively`,
		linalg::OperatorStoredMatrixMode	operator_stored_matrix_mode = `linalg::OperatorStoredMatrixMode::Off`
	)

inline

Member Function Documentation

◆ apply_impl()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::apply_impl	(	const SrcVectorType &	src,
		DstVectorType &	dst
	)

inline

Apply the operator: dst <- Op(src) (or additive, depending on mode).

Host-side responsibilities:

Handle replace/add mode initialization
Cache src/dst views for kernel capture
Build team policy
Dispatch to exactly one kernel variant based on kernel_path_
Optional halo communication accumulation

This is where the path decision now lives (host), so device code does not branch on kernel_path_ in the hot path.

◆ assemble_local_matrix()

template<typename ScalarT , int VecDim = 3>

dense::Mat< ScalarT, LocalMatrixDim, LocalMatrixDim > terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::assemble_local_matrix	(	const int	local_subdomain_id,
		const int	x_cell,
		const int	y_cell,
		const int	r_cell,
		const int	wedge
	)		const

inline

Assemble one wedge-local 18x18 matrix (slow path / on-demand assembly).

◆ assemble_trial_test_vecs()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::assemble_trial_test_vecs	(	const int	wedge,
		const dense::Vec< ScalarType, VecDim > &	quad_point,
		const ScalarType	quad_weight,
		const ScalarT	r_1,
		const ScalarT	r_2,
		dense::Vec< ScalarT, 3 >(*)	wedge_phy_surf[3],
		const dense::Vec< ScalarT, 6 > *	k_local_hex,
		const int	dimi,
		const int	dimj,
		dense::Mat< ScalarType, VecDim, VecDim > *	sym_grad_i,
		dense::Mat< ScalarType, VecDim, VecDim > *	sym_grad_j,
		ScalarType &	jdet_keval_quadweight
	)		const

inline

◆ column_grad_to_sym()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::column_grad_to_sym	(	const int	dim,
		const double	g0,
		const double	g1,
		const double	g2,
		double &	E00,
		double &	E11,
		double &	E22,
		double &	sym01,
		double &	sym02,
		double &	sym12,
		double &	gdd
	)		const

inline

Convert one gradient column into symmetric-gradient and div contributions.

Given grad(phi e_dim), this helper computes:

diagonal entries (E00,E11,E22)
off-diagonal symmetric entries (sym01,sym02,sym12)
divergence contribution gdd

The fast kernels use this repeatedly in fused operator application.

◆ get_domain()

template<typename ScalarT , int VecDim = 3>

const grid::shell::DistributedDomain & terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::get_domain ( ) const

inline

◆ get_grid()

template<typename ScalarT , int VecDim = 3>

grid::Grid3DDataVec< ScalarT, 3 > terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::get_grid ( )

inline

◆ get_local_matrix()

template<typename ScalarT , int VecDim = 3>

dense::Mat< ScalarT, LocalMatrixDim, LocalMatrixDim > terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::get_local_matrix	(	const int	local_subdomain_id,
		const int	x_cell,
		const int	y_cell,
		const int	r_cell,
		const int	wedge
	)		const

inline

◆ get_radii()

template<typename ScalarT , int VecDim = 3>

grid::Grid2DDataScalar< ScalarT > terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::get_radii ( ) const

inline

◆ get_stored_matrix_mode()

template<typename ScalarT , int VecDim = 3>

linalg::OperatorStoredMatrixMode terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::get_stored_matrix_mode ( )

inline

◆ has_flag()

template<typename ScalarT , int VecDim = 3>

bool terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::has_flag	(	const int	local_subdomain_id,
		const int	x_cell,
		const int	y_cell,
		const int	r_cell,
		grid::shell::ShellBoundaryFlag	flag
	)		const

inline

◆ k_grid_data()

template<typename ScalarT , int VecDim = 3>

const grid::Grid4DDataScalar< ScalarType > & terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::k_grid_data ( )

inline

◆ operator()()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::operator() ( const Team & team ) const

inline

Legacy generic team operator.

Kept for compatibility/debugging, but no longer used by apply_impl(). The host now dispatches directly to path-specific kernels.

This function still works, but it reintroduces a branch on kernel_path_ and should therefore be avoided in performance-critical use.

◆ set_boundary_conditions()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::set_boundary_conditions ( BoundaryConditions bcs )

inline

◆ set_diagonal()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::set_diagonal ( bool v )

inline

◆ set_local_matrix()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::set_local_matrix	(	const int	local_subdomain_id,
		const int	x_cell,
		const int	y_cell,
		const int	r_cell,
		const int	wedge,
		const dense::Mat< ScalarT, LocalMatrixDim, LocalMatrixDim > &	mat
	)		const

inline

◆ set_operator_apply_and_communication_modes()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::set_operator_apply_and_communication_modes	(	const linalg::OperatorApplyMode	operator_apply_mode,
		const linalg::OperatorCommunicationMode	operator_communication_mode
	)

inline

◆ set_stored_matrix_mode()

template<typename ScalarT , int VecDim = 3>

void terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::set_stored_matrix_mode	(	linalg::OperatorStoredMatrixMode	operator_stored_matrix_mode,
		int	level_range,
		grid::Grid4DDataScalar< ScalarType >	GCAElements
	)

inline

◆ team_shmem_size()

template<typename ScalarT , int VecDim = 3>

size_t terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::team_shmem_size ( const int ) const

inline

Team scratch memory size for fast paths.

Layout per team: [coords_sh | src_sh | k_sh | r_sh | padding]

Member Data Documentation

◆ LocalMatrixDim

template<typename ScalarT , int VecDim = 3>

constexpr int terra::fe::wedge::operators::shell::EpsilonDivDivKerngen< ScalarT, VecDim >::LocalMatrixDim = 18

staticconstexpr

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

epsilon_divdiv_kerngen.hpp

Public Types

Public Member Functions

Static Public Attributes

Detailed Description

IMPORTANT DESIGN CHANGE:

Member Typedef Documentation

◆ DstVectorType

◆ Grid4DDataLocalMatrices

◆ LocalMatrixStorage

◆ ScalarType

◆ SrcVectorType

◆ Team

Constructor & Destructor Documentation

◆ EpsilonDivDivKerngen()

Member Function Documentation

◆ apply_impl()

◆ assemble_local_matrix()

◆ assemble_trial_test_vecs()

◆ column_grad_to_sym()

◆ get_domain()

◆ get_grid()

◆ get_local_matrix()

◆ get_radii()

◆ get_stored_matrix_mode()

◆ has_flag()

◆ k_grid_data()

◆ operator()()

◆ set_boundary_conditions()

◆ set_diagonal()

◆ set_local_matrix()

◆ set_operator_apply_and_communication_modes()

◆ set_stored_matrix_mode()

◆ team_shmem_size()

Member Data Documentation

◆ LocalMatrixDim