Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

H(curl) basis and CEED_EVAL_CURL #1156

Merged
merged 8 commits into from
Apr 12, 2023
Merged

H(curl) basis and CEED_EVAL_CURL #1156

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
c4e3f59
Add H(curl) basis and CEED_EVAL_CURL operations
sebastiangrimberg Feb 15, 2023
15ad391
Update p-multigrid to work with H(div) and H(curl)
sebastiangrimberg Apr 3, 2023
352a5e7
Operator diagonal and full assembly support for H(curl) and H(div), o…
sebastiangrimberg Feb 24, 2023
5330800
Update function comments for CeedQFunctions related to H(div) and H(c…
sebastiangrimberg Feb 24, 2023
b5404d5
Add tests
sebastiangrimberg Feb 16, 2023
11b88dd
Update Julia/Python/Rust/Fortran bindings
sebastiangrimberg Apr 3, 2023
97c1c57
Address PR comments for high-level APIs
sebastiangrimberg Apr 6, 2023
b2e3f8e
Add Julia API tests
sebastiangrimberg Apr 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 5 additions & 39 deletions backends/blocked/ceed-blocked-operator.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo
case CEED_EVAL_INTERP:
case CEED_EVAL_GRAD:
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
Expand All @@ -92,8 +93,6 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo
CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
break;
case CEED_EVAL_CURL:
break; // Not implemented
}
}
return CEED_ERROR_SUCCESS;
Expand Down Expand Up @@ -227,33 +226,16 @@ static inline int CeedOperatorInputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFunc
CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * Q * size]));
break;
case CEED_EVAL_INTERP:
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i]));
break;
case CEED_EVAL_GRAD:
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i]));
break;
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i]));
break;
case CEED_EVAL_WEIGHT:
break; // No action
// LCOV_EXCL_START
case CEED_EVAL_CURL: {
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
Ceed ceed;
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
// LCOV_EXCL_STOP
}
}
}
return CEED_ERROR_SUCCESS;
Expand All @@ -280,36 +262,20 @@ static inline int CeedOperatorOutputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFun
case CEED_EVAL_NONE:
break; // No action
case CEED_EVAL_INTERP:
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(
CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i]));
break;
case CEED_EVAL_GRAD:
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(
CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i]));
break;
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(
CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
break;
// LCOV_EXCL_START
case CEED_EVAL_WEIGHT: {
Ceed ceed;
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
}
case CEED_EVAL_CURL: {
Ceed ceed;
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
// LCOV_EXCL_STOP
}
}
Expand Down
42 changes: 5 additions & 37 deletions backends/opt/ceed-opt-operator.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,7 @@ static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, bool i
case CEED_EVAL_INTERP:
case CEED_EVAL_GRAD:
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
Expand All @@ -95,8 +96,6 @@ static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, bool i
CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
break;
case CEED_EVAL_CURL:
break; // Not implemented
}
if (is_input && !!e_vecs[i]) {
CeedCallBackend(CeedVectorSetArray(e_vecs[i], CEED_MEM_HOST, CEED_COPY_VALUES, NULL));
Expand Down Expand Up @@ -248,39 +247,18 @@ static inline int CeedOperatorInputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunction
}
break;
case CEED_EVAL_INTERP:
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
if (!active_in) {
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp]));
}
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i]));
break;
case CEED_EVAL_GRAD:
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
if (!active_in) {
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp]));
}
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i]));
break;
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
if (!active_in) {
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp]));
}
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i]));
break;
case CEED_EVAL_WEIGHT:
break; // No action
// LCOV_EXCL_START
case CEED_EVAL_CURL: {
CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
Ceed ceed;
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
// LCOV_EXCL_STOP
}
}
}
return CEED_ERROR_SUCCESS;
Expand All @@ -306,16 +284,11 @@ static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunctio
case CEED_EVAL_NONE:
break; // No action
case CEED_EVAL_INTERP:
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i]));
break;
case CEED_EVAL_GRAD:
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i]));
break;
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i]));
CeedCallBackend(CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
break;
// LCOV_EXCL_START
case CEED_EVAL_WEIGHT: {
Expand All @@ -324,11 +297,6 @@ static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, CeedQFunctio
return CeedError(ceed, CEED_ERROR_BACKEND,
"CEED_EVAL_WEIGHT cannot be an output "
"evaluation mode");
}
case CEED_EVAL_CURL: {
Ceed ceed;
CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented");
// LCOV_EXCL_STOP
}
}
Expand Down
74 changes: 37 additions & 37 deletions backends/ref/ceed-ref-basis.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,12 @@
static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
Ceed ceed;
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
CeedInt dim, num_comp, num_nodes, num_qpts, Q_comp;
CeedInt dim, num_comp, q_comp, num_nodes, num_qpts;
CeedCallBackend(CeedBasisGetDimension(basis, &dim));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes));
CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, &Q_comp));
CeedTensorContract contract;
CeedCallBackend(CeedBasisGetTensorContract(basis, &contract));
const CeedInt add = (t_mode == CEED_TRANSPOSE);
Expand All @@ -46,8 +46,8 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo
}
bool tensor_basis;
CeedCallBackend(CeedBasisIsTensor(basis, &tensor_basis));
// Tensor basis
if (tensor_basis) {
// Tensor basis
CeedInt P_1d, Q_1d;
CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
Expand Down Expand Up @@ -191,36 +191,31 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo
}
} else {
// Non-tensor basis
CeedInt P = num_nodes, Q = num_qpts;
nbeams marked this conversation as resolved.
Show resolved Hide resolved
switch (eval_mode) {
// Interpolate to/from quadrature points
case CEED_EVAL_INTERP: {
CeedInt P = num_nodes, Q = Q_comp * num_qpts;
const CeedScalar *interp;
CeedCallBackend(CeedBasisGetInterp(basis, &interp));
if (t_mode == CEED_TRANSPOSE) {
P = Q_comp * num_qpts;
Q = num_nodes;
}
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, interp, t_mode, add, u, v));
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v));
} break;
// Evaluate the gradient to/from quadrature points
case CEED_EVAL_GRAD: {
CeedInt P = num_nodes, Q = num_qpts;
CeedInt dim_stride = num_qpts * num_comp * num_elem;
CeedInt grad_stride = num_qpts * num_nodes;
const CeedScalar *grad;
CeedCallBackend(CeedBasisGetGrad(basis, &grad));
if (t_mode == CEED_TRANSPOSE) {
P = num_qpts;
Q = num_nodes;
for (CeedInt d = 0; d < dim; d++) {
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, grad + d * grad_stride, t_mode, add, u + d * dim_stride, v));
}
} else {
for (CeedInt d = 0; d < dim; d++) {
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, grad + d * grad_stride, t_mode, add, u, v + d * dim_stride));
}
}
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v));
} break;
// Evaluate the divergence to/from the quadrature points
case CEED_EVAL_DIV: {
const CeedScalar *div;
CeedCallBackend(CeedBasisGetDiv(basis, &div));
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v));
} break;
// Evaluate the curl to/from the quadrature points
case CEED_EVAL_CURL: {
const CeedScalar *curl;
CeedCallBackend(CeedBasisGetCurl(basis, &curl));
CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v));
} break;
// Retrieve interpolation weights
case CEED_EVAL_WEIGHT: {
Expand All @@ -235,21 +230,7 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo
for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i];
}
} break;
// Evaluate the divergence to/from the quadrature points
case CEED_EVAL_DIV: {
CeedInt P = num_nodes, Q = num_qpts;
const CeedScalar *div;
CeedCallBackend(CeedBasisGetDiv(basis, &div));
if (t_mode == CEED_TRANSPOSE) {
P = num_qpts;
Q = num_nodes;
}
CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, div, t_mode, add, u, v));
} break;
// LCOV_EXCL_START
// Evaluate the curl to/from the quadrature points
case CEED_EVAL_CURL:
return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported");
// Take no action, BasisApply should not have been called
case CEED_EVAL_NONE:
return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
Expand Down Expand Up @@ -301,6 +282,25 @@ int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_node
return CEED_ERROR_SUCCESS;
}

//------------------------------------------------------------------------------
// Basis Create Non-Tensor H(curl)
//------------------------------------------------------------------------------
int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
Ceed ceed;
CeedCallBackend(CeedBasisGetCeed(basis, &ceed));

Ceed parent;
CeedCallBackend(CeedGetParent(ceed, &parent));
CeedTensorContract contract;
CeedCallBackend(CeedTensorContractCreate(parent, basis, &contract));
CeedCallBackend(CeedBasisSetTensorContract(basis, contract));

CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));

return CEED_ERROR_SUCCESS;
}

//------------------------------------------------------------------------------
// Basis Destroy Tensor
//------------------------------------------------------------------------------
Expand Down
Loading