Skip to content

Commit

Permalink
Better approach for libCEED QFunction organization with scalar-, vect…
Browse files Browse the repository at this point in the history 
…or-, and matrix-valued coefficient variants
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Sep 11, 2023
1 parent a786edc commit d020d9a
Show file tree
Hide file tree
Showing 19 changed files with 2,426 additions and 1,224 deletions.
56 changes: 39 additions & 17 deletions fem/ceed/integrators/convection/convection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,28 @@ namespace ceed
#ifdef MFEM_USE_CEED
struct ConvectionOperatorInfo : public OperatorInfo
{
ConvectionContext ctx;
ConvectionContext ctx = {0};
ConvectionOperatorInfo(const mfem::FiniteElementSpace &fes,
mfem::VectorCoefficient *VQ, double alpha,
bool use_bdr)
bool use_bdr = false, bool use_mf = false)
{
MFEM_VERIFY(VQ && VQ->GetVDim() == fes.GetMesh()->SpaceDimension(),
"Incorrect coefficient dimensions in ceed::ConvectionOperatorInfo!");
ctx.dim = fes.GetMesh()->Dimension() - use_bdr;
ctx.space_dim = fes.GetMesh()->SpaceDimension();
if (VectorConstantCoefficient *const_coeff =
dynamic_cast<VectorConstantCoefficient *>(VQ))
ctx.alpha = alpha;
if (!use_mf)
{
apply_func = ":f_apply_conv";
apply_qf = &f_apply_conv;
}
else
{
build_func = "";
build_qf = nullptr;
}
if (mfem::VectorConstantCoefficient *const_coeff =
dynamic_cast<mfem::VectorConstantCoefficient *>(VQ))
{
const int vdim = VQ->GetVDim();
MFEM_VERIFY(vdim <= LIBCEED_CONV_COEFF_COMP_MAX,
Expand All @@ -45,20 +56,31 @@ struct ConvectionOperatorInfo : public OperatorInfo
{
ctx.coeff[i] = val[i];
}
if (!use_mf)
{
build_func = ":f_build_conv_const";
build_qf = &f_build_conv_const;
}
else
{
apply_func = ":f_apply_conv_mf_const";
apply_qf = &f_apply_conv_mf_const;
}
}
else
{
if (!use_mf)
{
build_func = ":f_build_conv_quad";
build_qf = &f_build_conv_quad;
}
else
{
apply_func = ":f_apply_conv_mf_quad";
apply_qf = &f_apply_conv_mf_quad;
}
}
ctx.alpha = alpha;

header = "/integrators/convection/convection_qf.h";
build_func_const = ":f_build_conv_const";
build_qf_const = &f_build_conv_const;
build_func_quad = ":f_build_conv_quad";
build_qf_quad = &f_build_conv_quad;
apply_func = ":f_apply_conv";
apply_qf = &f_apply_conv;
apply_func_mf_const = ":f_apply_conv_mf_const";
apply_qf_mf_const = &f_apply_conv_mf_const;
apply_func_mf_quad = ":f_apply_conv_mf_quad";
apply_qf_mf_quad = &f_apply_conv_mf_quad;
trial_op = EvalMode::Grad;
test_op = EvalMode::Interp;
qdatasize = ctx.dim;
Expand Down Expand Up @@ -89,7 +111,7 @@ MFConvectionIntegrator::MFConvectionIntegrator(
const bool use_bdr)
{
#ifdef MFEM_USE_CEED
ConvectionOperatorInfo info(fes, VQ, alpha, use_bdr);
ConvectionOperatorInfo info(fes, VQ, alpha, use_bdr, true);
Assemble(integ, info, fes, VQ, use_bdr, true);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
Expand Down
21 changes: 11 additions & 10 deletions fem/ceed/integrators/convection/convection_qf.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,11 @@

#define LIBCEED_CONV_COEFF_COMP_MAX 3

/// A structure used to pass additional data to f_build_conv and f_apply_conv
struct ConvectionContext
{
CeedInt dim, space_dim;
CeedScalar coeff[LIBCEED_CONV_COEFF_COMP_MAX];
CeedScalar alpha;
CeedScalar coeff[LIBCEED_CONV_COEFF_COMP_MAX];
};

/// libCEED QFunction for building quadrature data for a convection operator
Expand All @@ -34,7 +33,7 @@ CEED_QFUNCTION(f_build_conv_const)(void *ctx, CeedInt Q,
// in[0] is Jacobians with shape [dim, ncomp=space_dim, Q]
// in[1] is quadrature weights, size (Q)
//
// At every quadrature point, compute and store qw * α * c^T adj(J)^T.
// At every quadrature point, compute and store qw * α * c^T adj(J)^T
const CeedScalar alpha = bc->alpha;
const CeedScalar *coeff = bc->coeff;
const CeedScalar *J = in[0], *qw = in[1];
Expand Down Expand Up @@ -77,7 +76,7 @@ CEED_QFUNCTION(f_build_conv_const)(void *ctx, CeedInt Q,
}

/// libCEED QFunction for building quadrature data for a convection operator
/// with a coefficient evaluated at quadrature points.
/// with a coefficient evaluated at quadrature points
CEED_QFUNCTION(f_build_conv_quad)(void *ctx, CeedInt Q,
const CeedScalar *const *in,
CeedScalar *const *out)
Expand All @@ -87,7 +86,7 @@ CEED_QFUNCTION(f_build_conv_quad)(void *ctx, CeedInt Q,
// in[1] is Jacobians with shape [dim, ncomp=space_dim, Q]
// in[2] is quadrature weights, size (Q)
//
// At every quadrature point, compute and store qw * α * c^T adj(J)^T.
// At every quadrature point, compute and store qw * α * c^T adj(J)^T
const CeedScalar alpha = bc->alpha;
const CeedScalar *c = in[0], *J = in[1], *qw = in[2];
CeedScalar *qd = out[0];
Expand Down Expand Up @@ -127,7 +126,7 @@ CEED_QFUNCTION(f_build_conv_quad)(void *ctx, CeedInt Q,
return 0;
}

/// libCEED QFunction for applying a conv operator
/// libCEED QFunction for applying a convection operator
CEED_QFUNCTION(f_apply_conv)(void *ctx, CeedInt Q,
const CeedScalar *const *in,
CeedScalar *const *out)
Expand Down Expand Up @@ -166,7 +165,8 @@ CEED_QFUNCTION(f_apply_conv)(void *ctx, CeedInt Q,
return 0;
}

/// libCEED QFunction for applying a conv operator
/// libCEED QFunction for applying a convection operator with a constant
/// coefficient
CEED_QFUNCTION(f_apply_conv_mf_const)(void *ctx, CeedInt Q,
const CeedScalar *const *in,
CeedScalar *const *out)
Expand All @@ -177,7 +177,7 @@ CEED_QFUNCTION(f_apply_conv_mf_const)(void *ctx, CeedInt Q,
// in[2] is quadrature weights, size (Q)
// out[0] has shape [ncomp=1, Q]
//
// At every quadrature point, compute qw * α * c^T adj(J)^T.
// At every quadrature point, compute qw * α * c^T adj(J)^T
const CeedScalar alpha = bc->alpha;
const CeedScalar *coeff = bc->coeff;
const CeedScalar *ug = in[0], *J = in[1], *qw = in[2];
Expand Down Expand Up @@ -235,7 +235,8 @@ CEED_QFUNCTION(f_apply_conv_mf_const)(void *ctx, CeedInt Q,
return 0;
}

/// libCEED QFunction for applying a conv operator
/// libCEED QFunction for applying a convection operator with a coefficient
/// evaluated at quadrature points
CEED_QFUNCTION(f_apply_conv_mf_quad)(void *ctx, CeedInt Q,
const CeedScalar *const *in,
CeedScalar *const *out)
Expand All @@ -247,7 +248,7 @@ CEED_QFUNCTION(f_apply_conv_mf_quad)(void *ctx, CeedInt Q,
// in[3] is quadrature weights, size (Q)
// out[0] has shape [ncomp=1, Q]
//
// At every quadrature point, compute qw * α * c^T adj(J)^T.
// At every quadrature point, compute qw * α * c^T adj(J)^T
const CeedScalar alpha = bc->alpha;
const CeedScalar *ug = in[0], *c = in[1], *J = in[2], *qw = in[3];
CeedScalar *vg = out[0];
Expand Down
138 changes: 106 additions & 32 deletions fem/ceed/integrators/curlcurl/curlcurl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,10 @@ namespace ceed
#ifdef MFEM_USE_CEED
struct CurlCurlOperatorInfo : public OperatorInfo
{
CurlCurlContext ctx;
CurlCurlContext ctx = {0};
template <typename CoeffType>
CurlCurlOperatorInfo(const mfem::FiniteElementSpace &fes, CoeffType *Q,
bool use_bdr)
bool use_bdr = false, bool use_mf = false)
{
MFEM_VERIFY(fes.GetVDim() == 1,
"libCEED interface for vector FE does not support VDim > 1!");
Expand All @@ -37,65 +37,116 @@ struct CurlCurlOperatorInfo : public OperatorInfo
"CurlCurlIntegrator requires dim == 2 or dim == 3!");
ctx.space_dim = fes.GetMesh()->SpaceDimension();
ctx.curl_dim = (ctx.dim < 3) ? 1 : ctx.dim;
if (!use_mf)
{
apply_func = ":f_apply_curlcurl";
apply_qf = &f_apply_curlcurl;
}
else
{
build_func = "";
build_qf = nullptr;
}
if (Q == nullptr)
{
ctx.coeff_comp = 1;
ctx.coeff[0] = 1.0;
if (!use_mf)
{
build_func = ":f_build_curlcurl_const_scalar";
build_qf = &f_build_curlcurl_const_scalar;
}
else
{
apply_func = ":f_apply_curlcurl_mf_const_scalar";
apply_qf = &f_apply_curlcurl_mf_const_scalar;
}
}
else
{
InitCoefficient(*Q);
InitCoefficient(*Q, use_mf);
}

header = "/integrators/curlcurl/curlcurl_qf.h";
build_func_const = ":f_build_curlcurl_const";
build_qf_const = &f_build_curlcurl_const;
build_func_quad = ":f_build_curlcurl_quad";
build_qf_quad = &f_build_curlcurl_quad;
apply_func = ":f_apply_curlcurl";
apply_qf = &f_apply_curlcurl;
apply_func_mf_const = ":f_apply_curlcurl_mf_const";
apply_qf_mf_const = &f_apply_curlcurl_mf_const;
apply_func_mf_quad = ":f_apply_curlcurl_mf_quad";
apply_qf_mf_quad = &f_apply_curlcurl_mf_quad;
trial_op = EvalMode::Curl;
test_op = EvalMode::Curl;
qdatasize = (ctx.curl_dim * (ctx.curl_dim + 1)) / 2;
}
void InitCoefficient(mfem::Coefficient &Q)
void InitCoefficient(mfem::Coefficient &Q, bool use_mf)
{
ctx.coeff_comp = 1;
if (ConstantCoefficient *const_coeff =
dynamic_cast<ConstantCoefficient *>(&Q))
if (mfem::ConstantCoefficient *const_coeff =
dynamic_cast<mfem::ConstantCoefficient *>(&Q))
{
ctx.coeff[0] = const_coeff->constant;
if (!use_mf)
{
build_func = ":f_build_curlcurl_const_scalar";
build_qf = &f_build_curlcurl_const_scalar;
}
else
{
apply_func = ":f_apply_curlcurl_mf_const_scalar";
apply_qf = &f_apply_curlcurl_mf_const_scalar;
}
}
else
{
if (!use_mf)
{
build_func = ":f_build_curlcurl_quad_scalar";
build_qf = &f_build_curlcurl_quad_scalar;
}
else
{
apply_func = ":f_apply_curlcurl_mf_quad_scalar";
apply_qf = &f_apply_curlcurl_mf_quad_scalar;
}
}
}
void InitCoefficient(mfem::VectorCoefficient &VQ)
void InitCoefficient(mfem::VectorCoefficient &VQ, bool use_mf)
{
const int vdim = VQ.GetVDim();
ctx.coeff_comp = vdim;
if (VectorConstantCoefficient *const_coeff =
dynamic_cast<VectorConstantCoefficient *>(&VQ))
if (mfem::VectorConstantCoefficient *const_coeff =
dynamic_cast<mfem::VectorConstantCoefficient *>(&VQ))
{
MFEM_VERIFY(ctx.coeff_comp <= LIBCEED_CURLCURL_COEFF_COMP_MAX,
const int vdim = VQ.GetVDim();
MFEM_VERIFY(vdim <= LIBCEED_CURLCURL_COEFF_COMP_MAX,
"VectorCoefficient dimension exceeds context storage!");
const mfem::Vector &val = const_coeff->GetVec();
for (int i = 0; i < vdim; i++)
{
ctx.coeff[i] = val[i];
}
if (!use_mf)
{
build_func = ":f_build_curlcurl_const_vector";
build_qf = &f_build_curlcurl_const_vector;
}
else
{
apply_func = ":f_apply_curlcurl_mf_const_vector";
apply_qf = &f_apply_curlcurl_mf_const_vector;
}
}
else
{
if (!use_mf)
{
build_func = ":f_build_curlcurl_quad_vector";
build_qf = &f_build_curlcurl_quad_vector;
}
else
{
apply_func = ":f_apply_curlcurl_mf_quad_vector";
apply_qf = &f_apply_curlcurl_mf_quad_vector;
}
}
}
void InitCoefficient(mfem::MatrixCoefficient &MQ)
void InitCoefficient(mfem::MatrixCoefficient &MQ, bool use_mf)
{
// Assumes matrix coefficient is symmetric
const int vdim = MQ.GetVDim();
ctx.coeff_comp = (vdim * (vdim + 1)) / 2;
if (MatrixConstantCoefficient *const_coeff =
dynamic_cast<MatrixConstantCoefficient *>(&MQ))
if (mfem::MatrixConstantCoefficient *const_coeff =
dynamic_cast<mfem::MatrixConstantCoefficient *>(&MQ))
{
MFEM_VERIFY(ctx.coeff_comp <= LIBCEED_CURLCURL_COEFF_COMP_MAX,
const int vdim = MQ.GetVDim();
MFEM_VERIFY((vdim * (vdim + 1)) / 2 <= LIBCEED_CURLCURL_COEFF_COMP_MAX,
"MatrixCoefficient dimensions exceed context storage!");
const mfem::DenseMatrix &val = const_coeff->GetMatrix();
for (int j = 0; j < vdim; j++)
Expand All @@ -106,6 +157,29 @@ struct CurlCurlOperatorInfo : public OperatorInfo
ctx.coeff[idx] = val(i, j);
}
}
if (!use_mf)
{
build_func = ":f_build_curlcurl_const_matrix";
build_qf = &f_build_curlcurl_const_matrix;
}
else
{
apply_func = ":f_apply_curlcurl_mf_const_matrix";
apply_qf = &f_apply_curlcurl_mf_const_matrix;
}
}
else
{
if (!use_mf)
{
build_func = ":f_build_curlcurl_quad_matrix";
build_qf = &f_build_curlcurl_quad_matrix;
}
else
{
apply_func = ":f_apply_curlcurl_mf_quad_matrix";
apply_qf = &f_apply_curlcurl_mf_quad_matrix;
}
}
}
};
Expand Down Expand Up @@ -134,7 +208,7 @@ MFCurlCurlIntegrator::MFCurlCurlIntegrator(
const bool use_bdr)
{
#ifdef MFEM_USE_CEED
CurlCurlOperatorInfo info(fes, Q, use_bdr);
CurlCurlOperatorInfo info(fes, Q, use_bdr, true);
Assemble(integ, info, fes, Q, use_bdr, true);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
Expand Down

0 comments on commit d020d9a

Please sign in to comment.