From 052a4caf32ba065d5f89f3a216eb78801e55b245 Mon Sep 17 00:00:00 2001 From: "Kenneth E. Jansen" Date: Mon, 24 Jul 2023 23:56:12 +0000 Subject: [PATCH 01/15] breaking lots of DL stuff to test PetscFE --- examples/fluids/navierstokes.h | 2 +- examples/fluids/problems/sgs_dd_model.c | 4 +- examples/fluids/src/differential_filter.c | 4 +- examples/fluids/src/grid_anisotropy_tensor.c | 7 +- examples/fluids/src/setupdm.c | 35 ++- examples/fluids/src/setuplibceed.c | 238 +++++++++++++++++- .../fluids/src/strong_boundary_conditions.c | 5 +- examples/fluids/src/turb_spanstats.c | 4 +- .../fluids/src/velocity_gradient_projection.c | 4 +- 9 files changed, 276 insertions(+), 27 deletions(-) diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index 30d96eeeed..458bb2fb0d 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -322,7 +322,7 @@ PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLab CeedElemRestriction *elem_restr); // Utility function to get Ceed Restriction for each domain -PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, +PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, CeedInt Q_dim, CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, CeedElemRestriction *elem_restr_qd_i); diff --git a/examples/fluids/problems/sgs_dd_model.c b/examples/fluids/problems/sgs_dd_model.c index 33eac19293..1e76b3dad4 100644 --- a/examples/fluids/problems/sgs_dd_model.c +++ b/examples/fluids/problems/sgs_dd_model.c @@ -86,8 +86,8 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_grad_velo); CeedElemRestrictionGetNumComponents(elem_restr_grad_velo, &num_comp_grad_velo); } - - PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, num_qpts_1d, 0, &elem_restr_sgs, NULL, NULL)); +// 0 added FIXME + PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, num_qpts_1d, 0,0, &elem_restr_sgs, NULL, NULL)); CeedElemRestrictionCreateVector(elem_restr_sgs, &sgs_dd_data->sgs_nodal_ceed, NULL); // -- Create inverse multiplicity for correcting nodal assembly diff --git a/examples/fluids/src/differential_filter.c b/examples/fluids/src/differential_filter.c index 7a05ff5d19..b0970c487c 100644 --- a/examples/fluids/src/differential_filter.c +++ b/examples/fluids/src/differential_filter.c @@ -64,8 +64,8 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData char field_name[PETSC_MAX_PATH_LEN]; CeedElemRestriction elem_restr_filter; CeedBasis basis_filter; - - PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, num_qpts_1d, 0, &elem_restr_filter, NULL, NULL)); +// 0 added FIXME + PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, num_qpts_1d, 0,0, &elem_restr_filter, NULL, NULL)); CeedBasisCreateTensorH1Lagrange(ceed, dim, diff_filter->num_field_components[i], num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_filter); PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i)); diff --git a/examples/fluids/src/grid_anisotropy_tensor.c b/examples/fluids/src/grid_anisotropy_tensor.c index 6cceabc156..f8c250a4b5 100644 --- a/examples/fluids/src/grid_anisotropy_tensor.c +++ b/examples/fluids/src/grid_anisotropy_tensor.c @@ -59,8 +59,9 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); +//HACK to get to compile but won't work until code changed to get basis from PetscFE added extra 0 to arg list FIXME PetscCall( - GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, num_qpts_1d, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); + GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, num_qpts_1d, 0,grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); CeedBasisCreateTensorH1Lagrange(ceed, dim, grid_aniso_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grid_aniso); // -- Build RHS operator @@ -145,8 +146,8 @@ PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User use *num_comp_aniso = 7; CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); - - PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_qpts_1d, *num_comp_aniso, NULL, NULL, elem_restr_grid_aniso)); +// here too FIXME + PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_qpts_1d, 0,*num_comp_aniso, NULL, NULL, elem_restr_grid_aniso)); // -- Build collocation operator CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc); diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index aa9912fb5d..5369626721 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -41,14 +41,40 @@ PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, MatType mat_type, V PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc, Physics phys) { PetscFunctionBeginUser; { + PetscBool is_simplex = PETSC_TRUE; + + // Check if simplex or tensor-product mesh + PetscCall(DMPlexIsSimplex(dm, &is_simplex)); // Configure the finite element space and boundary conditions PetscFE fe; PetscInt num_comp_q = 5; DMLabel label; - PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, is_simplex, degree, PETSC_DECIDE, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Q")); PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); PetscCall(DMCreateDS(dm)); + // Project coordinates to enrich quadrature space + { + DM dm_coord; + PetscDS ds_coord; + PetscFE fe_coord_current, fe_coord_new; + PetscDualSpace fe_coord_dual_space; + PetscInt fe_coord_order, num_comp_coord; + + PetscCall(DMGetCoordinateDM(dm, &dm_coord)); + PetscCall(DMGetCoordinateDim(dm, &num_comp_coord)); + PetscCall(DMGetRegionDS(dm_coord, NULL, NULL, &ds_coord,NULL)); + PetscCall(PetscDSGetDiscretization(ds_coord, 0, (PetscObject *)&fe_coord_current)); + PetscCall(PetscFEGetDualSpace(fe_coord_current, &fe_coord_dual_space)); + PetscCall(PetscDualSpaceGetOrder(fe_coord_dual_space, &fe_coord_order)); + + // Create FE for coordinates + fe_coord_order = 1 ; + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_coord, is_simplex, fe_coord_order, degree, &fe_coord_new)); + PetscCall(DMProjectCoordinates(dm, fe_coord_new)); + PetscCall(PetscFEDestroy(&fe_coord_new)); + } + PetscCall(DMGetLabel(dm, "Face Sets", &label)); PetscCall(DMPlexLabelComplete(dm, label)); // Set wall BCs @@ -77,6 +103,13 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc } PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); + if (!is_simplex) { + DM dm_coord; + PetscCall(DMGetCoordinateDM(dm, &dm_coord)); + PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); + PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); + } + PetscCall(PetscFEDestroy(&fe)); } diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index 0871d55234..ea028b05cb 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -12,6 +12,8 @@ #include #include "../navierstokes.h" +#include + // Utility function - essential BC dofs are encoded in closure indices as -(i+1). PetscInt Involute(PetscInt i) { return i >= 0 ? i : -(i + 1); } @@ -35,17 +37,18 @@ PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLab // Utility function to get Ceed Restriction for each domain PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, - CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, + CeedInt Q_dim, CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, CeedElemRestriction *elem_restr_qd_i) { DM dm_coord; - CeedInt Q_dim, loc_num_elem; +// now passed in/out CeedInt Q_dim + CeedInt loc_num_elem; PetscInt dim; CeedElemRestriction elem_restr_tmp; PetscFunctionBeginUser; PetscCall(DMGetDimension(dm, &dim)); dim -= height; - Q_dim = CeedIntPow(Q, dim); +// now passed in/out for simplex set by PetscFE Q_dim = CeedIntPow(Q, dim); PetscCall(CreateRestrictionFromPlex(ceed, dm, height, domain_label, label_value, dm_field, &elem_restr_tmp)); if (elem_restr_q) *elem_restr_q = elem_restr_tmp; if (elem_restr_x) { @@ -78,11 +81,11 @@ PetscErrorCode AddBCSubOperator(Ceed ceed, DM dm, CeedData ceed_data, DMLabel do CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); // ---- CEED Restriction - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, Q_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, Q_sur, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, &elem_restr_qd_i_sur)); if (jac_data_size_sur > 0) { // State-dependent data will be passed from residual to Jacobian. This will be collocated. - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, Q_sur, jac_data_size_sur, NULL, NULL, &elem_restr_jd_i_sur)); + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, Q_sur, num_qpts_sur, jac_data_size_sur, NULL, NULL, &elem_restr_jd_i_sur)); CeedElemRestrictionCreateVector(elem_restr_jd_i_sur, &jac_data_sur, NULL); } else { elem_restr_jd_i_sur = NULL; @@ -213,6 +216,205 @@ PetscErrorCode SetupBCQFunctions(Ceed ceed, PetscInt dim_sur, PetscInt num_comp_ PetscFunctionReturn(PETSC_SUCCESS); } + +// ----------------------------------------------------------------------------- +// Convert DM field to DS field +// ----------------------------------------------------------------------------- +PetscErrorCode DMFieldToDSField(DM dm, DMLabel domain_label, PetscInt dm_field, PetscInt *ds_field) { + PetscDS ds; + IS field_is; + const PetscInt *fields; + PetscInt num_fields; + + PetscFunctionBeginUser; + + // Translate dm_field to ds_field + PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL)); + PetscCall(ISGetIndices(field_is, &fields)); + PetscCall(ISGetSize(field_is, &num_fields)); + for (PetscInt i = 0; i < num_fields; i++) { + if (dm_field == fields[i]) { + *ds_field = i; + break; + } + } + PetscCall(ISRestoreIndices(field_is, &fields)); + + if (*ds_field == -1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field); + + PetscFunctionReturn(0); +} + + +// ----------------------------------------------------------------------------- +// Utility function - convert from DMPolytopeType to CeedElemTopology +// ----------------------------------------------------------------------------- +CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { + switch (cell_type) { + case DM_POLYTOPE_TRIANGLE: + return CEED_TOPOLOGY_TRIANGLE; + case DM_POLYTOPE_QUADRILATERAL: + return CEED_TOPOLOGY_QUAD; + case DM_POLYTOPE_TETRAHEDRON: + return CEED_TOPOLOGY_TET; + case DM_POLYTOPE_HEXAHEDRON: + return CEED_TOPOLOGY_HEX; + default: + return 0; + } +} + + +// ----------------------------------------------------------------------------- +// Create libCEED Basis from PetscTabulation +// ----------------------------------------------------------------------------- +PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt face, PetscFE fe, + PetscTabulation basis_tabulation, PetscQuadrature quadrature, CeedBasis *basis) { + PetscInt first_point; + PetscInt ids[1] = {label_value}; + DMLabel depth_label; + DMPolytopeType cell_type; + CeedElemTopology elem_topo; + PetscScalar *q_points, *interp, *grad; + const PetscScalar *q_weights; + PetscDualSpace dual_space; + PetscInt num_dual_basis_vectors; + PetscInt dim, num_comp, P, Q; + + PetscFunctionBeginUser; + + // General basis information + PetscCall(PetscFEGetSpatialDimension(fe, &dim)); + PetscCall(PetscFEGetNumComponents(fe, &num_comp)); + PetscCall(PetscFEGetDualSpace(fe, &dual_space)); + PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); + P = num_dual_basis_vectors / num_comp; + + // Use depth label if no domain label present + if (!domain_label) { + PetscInt depth; + + PetscCall(DMPlexGetDepth(dm, &depth)); + PetscCall(DMPlexGetDepthLabel(dm, &depth_label)); + ids[0] = depth - height; + } + + // Get cell interp, grad, and quadrature data + PetscCall(DMGetFirstLabeledPoint(dm, dm, domain_label ? domain_label : depth_label, 1, ids, height, &first_point, NULL)); + PetscCall(DMPlexGetCellType(dm, first_point, &cell_type)); + elem_topo = ElemTopologyP2C(cell_type); + if (!elem_topo) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMPlex topology not supported"); + { + size_t q_points_size; + const PetscScalar *q_points_petsc; + PetscInt q_dim; + + PetscCall(PetscQuadratureGetData(quadrature, &q_dim, NULL, &Q, &q_points_petsc, &q_weights)); + q_points_size = Q * dim * sizeof(CeedScalar); + PetscCall(PetscCalloc(q_points_size, &q_points)); + for (PetscInt q = 0; q < Q; q++) { + for (PetscInt d = 0; d < q_dim; d++) q_points[q * dim + d] = q_points_petsc[q * q_dim + d]; + } + } + + // Convert to libCEED orientation + { + PetscBool is_simplex = PETSC_FALSE; + IS permutation = NULL; + const PetscInt *permutation_indices; + + PetscCall(DMPlexIsSimplex(dm, &is_simplex)); + if (!is_simplex) { + PetscSection section; + + // -- Get permutation + PetscCall(DMGetLocalSection(dm, §ion)); + PetscCall(PetscSectionGetClosurePermutation(section, (PetscObject)dm, dim, num_comp * P, &permutation)); + PetscCall(ISGetIndices(permutation, &permutation_indices)); + } + + // -- Copy interp, grad matrices + PetscCall(PetscCalloc(P * Q * sizeof(CeedScalar), &interp)); + PetscCall(PetscCalloc(P * Q * dim * sizeof(CeedScalar), &grad)); + const CeedInt c = 0; + for (CeedInt q = 0; q < Q; q++) { + for (CeedInt p_ceed = 0; p_ceed < P; p_ceed++) { + CeedInt p_petsc = is_simplex ? (p_ceed * num_comp) : permutation_indices[p_ceed * num_comp]; + + interp[q * P + p_ceed] = basis_tabulation->T[0][((face * Q + q) * P * num_comp + p_petsc) * num_comp + c]; + for (CeedInt d = 0; d < dim; d++) { + grad[(d * Q + q) * P + p_ceed] = basis_tabulation->T[1][(((face * Q + q) * P * num_comp + p_petsc) * num_comp + c) * dim + d]; + } + } + } + + // -- Cleanup + if (permutation) PetscCall(ISRestoreIndices(permutation, &permutation_indices)); + PetscCall(ISDestroy(&permutation)); + } + + // Finally, create libCEED basis + CeedBasisCreateH1(ceed, elem_topo, num_comp, P, Q, interp, grad, q_points, q_weights, basis); + PetscCall(PetscFree(q_points)); + PetscCall(PetscFree(interp)); + PetscCall(PetscFree(grad)); + + PetscFunctionReturn(0); +} + + +// ----------------------------------------------------------------------------- +// Get CEED Basis from DMPlex +// ----------------------------------------------------------------------------- +PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedQuadMode quadMode, + CeedBasis *basis) { + PetscDS ds; + PetscFE fe; + PetscQuadrature quadrature; + PetscBool is_simplex = PETSC_TRUE; + PetscInt ds_field = -1; + + PetscFunctionBeginUser; + + // Get element information + PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL)); + PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field)); + PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe)); + PetscCall(PetscFEGetHeightSubspace(fe, height, &fe)); + PetscCall(PetscFEGetQuadrature(fe, &quadrature)); + + // Check if simplex or tensor-product mesh + PetscCall(DMPlexIsSimplex(dm, &is_simplex)); + + // Build libCEED basis + if (is_simplex) { + PetscTabulation basis_tabulation; + PetscInt num_derivatives = 1, face = 0; + + PetscCall(PetscFEGetCellTabulation(fe, num_derivatives, &basis_tabulation)); + PetscCall(BasisCreateFromTabulation(ceed, dm, domain_label, label_value, height, face, fe, basis_tabulation, quadrature, basis)); + } else { + PetscDualSpace dual_space; + PetscInt num_dual_basis_vectors; + PetscInt dim, num_comp, P, Q; + + PetscCall(PetscFEGetSpatialDimension(fe, &dim)); + PetscCall(PetscFEGetNumComponents(fe, &num_comp)); + PetscCall(PetscFEGetDualSpace(fe, &dual_space)); + PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); + P = num_dual_basis_vectors / num_comp; + PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &Q, NULL, NULL)); + + CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim)); + CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim)); + + CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, quadMode, basis); + } + + PetscFunctionReturn(0); +} + + PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { PetscFunctionBeginUser; @@ -224,22 +426,28 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App P = app_ctx->degree + 1, Q = P + app_ctx->q_extra; CeedElemRestriction elem_restr_jd_i; CeedVector jac_data; + CeedInt num_qpts; // ----------------------------------------------------------------------------- // CEED Bases // ----------------------------------------------------------------------------- - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_q, P, Q, CEED_GAUSS, &ceed_data->basis_q); - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, &ceed_data->basis_x); - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, CEED_GAUSS_LOBATTO, &ceed_data->basis_xc); + DM dm_coord; + PetscCall(DMGetCoordinateDM(dm, &dm_coord)); + + PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, CEED_GAUSS, &ceed_data->basis_q)); + PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, CEED_GAUSS, &ceed_data->basis_x)); + PetscCall(CeedBasisCreateProjection(ceed_data->basis_x, ceed_data->basis_q, &ceed_data->basis_xc)); + // --- Get number of quadrature points for the boundaries + CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts); // ----------------------------------------------------------------------------- // CEED Restrictions // ----------------------------------------------------------------------------- // -- Create restriction - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, q_data_size_vol, &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, num_qpts, q_data_size_vol, &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, &ceed_data->elem_restr_qd_i)); - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, jac_data_size_vol, NULL, NULL, &elem_restr_jd_i)); + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, num_qpts, jac_data_size_vol, NULL, NULL, &elem_restr_jd_i)); // -- Create E vectors CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, NULL); @@ -403,9 +611,13 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App // ----------------------------------------------------------------------------- // CEED Bases // ----------------------------------------------------------------------------- - CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_q, P_sur, Q_sur, CEED_GAUSS, &ceed_data->basis_q_sur); - CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, Q_sur, CEED_GAUSS, &ceed_data->basis_x_sur); - CeedBasisCreateTensorH1Lagrange(ceed, dim_sur, num_comp_x, 2, P_sur, CEED_GAUSS_LOBATTO, &ceed_data->basis_xc_sur); + + DMLabel label = 0; + PetscInt face_id = 0; + PetscInt field = 0; // Still want the normal, default field + PetscCall(CreateBasisFromPlex(ceed, dm, label, face_id, height, field, CEED_GAUSS, &ceed_data->basis_q_sur)); + PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, face_id, height, field, CEED_GAUSS, &ceed_data->basis_x_sur)); + PetscCall(CeedBasisCreateProjection(ceed_data->basis_x_sur, ceed_data->basis_q_sur, &ceed_data->basis_xc_sur)); // ----------------------------------------------------------------------------- // CEED QFunctions diff --git a/examples/fluids/src/strong_boundary_conditions.c b/examples/fluids/src/strong_boundary_conditions.c index 26b136ece9..10aabc0e97 100644 --- a/examples/fluids/src/strong_boundary_conditions.c +++ b/examples/fluids/src/strong_boundary_conditions.c @@ -29,6 +29,9 @@ PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, CeedData ceed_data, DM dm, Problem // Basis CeedInt height = 1; PetscCall(CeedBasisCreateProjection(ceed_data->basis_x_sur, ceed_data->basis_q_sur, &basis_x_to_q_sur)); + // --- Get number of quadrature points for the boundaries + CeedInt num_qpts_sur; + CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); // Setup QFunction CeedQFunctionCreateInterior(ceed, 1, SetupStrongBC, SetupStrongBC_loc, &qf_setup); @@ -43,7 +46,7 @@ PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, CeedData ceed_data, DM dm, Problem // Compute contribution on each boundary face for (CeedInt i = 0; i < bc->num_inflow; i++) { // -- Restrictions - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, Q_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, Q_sur, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, &elem_restr_qd_sur)); CeedElemRestrictionCreateVector(elem_restr_q_sur, &multiplicity, NULL); CeedElemRestrictionGetMultiplicity(elem_restr_q_sur, multiplicity); diff --git a/examples/fluids/src/turb_spanstats.c b/examples/fluids/src/turb_spanstats.c index 934d458d09..86150c2086 100644 --- a/examples/fluids/src/turb_spanstats.c +++ b/examples/fluids/src/turb_spanstats.c @@ -187,8 +187,8 @@ PetscErrorCode SpanStatsSetupDataCreate(Ceed ceed, User user, CeedData ceed_data PetscCall(DMGetDimension(dm, &dim)); CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &Q); CeedBasisGetNumNodes1D(ceed_data->basis_q, &P); - - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, problem->q_data_size_sur, &(*stats_data)->elem_restr_parent_stats, +// inserted this zero just to get it to compile FIXME + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, 0, problem->q_data_size_sur, &(*stats_data)->elem_restr_parent_stats, &(*stats_data)->elem_restr_parent_x, &(*stats_data)->elem_restr_parent_qd)); CeedElemRestrictionGetNumComponents((*stats_data)->elem_restr_parent_x, &num_comp_x); CeedElemRestrictionCreateVector((*stats_data)->elem_restr_parent_x, &(*stats_data)->x_coord, NULL); diff --git a/examples/fluids/src/velocity_gradient_projection.c b/examples/fluids/src/velocity_gradient_projection.c index 7a84fd2763..2b038bc6b8 100644 --- a/examples/fluids/src/velocity_gradient_projection.c +++ b/examples/fluids/src/velocity_gradient_projection.c @@ -69,8 +69,8 @@ PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, User user, CeedData ce CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); - - PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, num_qpts_1d, q_data_size, &elem_restr_grad_velo, NULL, NULL)); +// inserted zero just to compile FIXME + PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, num_qpts_1d, 0, q_data_size, &elem_restr_grad_velo, NULL, NULL)); CeedBasisCreateTensorH1Lagrange(ceed, dim, grad_velo_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grad_velo); From fbcaac505da2abe2059dd662cecfecc1a8686212 Mon Sep 17 00:00:00 2001 From: "Kenneth E. Jansen" Date: Tue, 25 Jul 2023 14:33:51 +0000 Subject: [PATCH 02/15] fixed missing header --- examples/fluids/src/setupdm.c | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index 5369626721..ba200a3a4a 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -10,6 +10,7 @@ #include #include +#include #include "../navierstokes.h" #include "../problems/stg_shur14.h" From 44dc348c69e5b78680b3a7345e44f5e97f381bd9 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 25 Jul 2023 10:44:24 -0600 Subject: [PATCH 03/15] fluids: Remove unused `Involute` function --- examples/fluids/navierstokes.h | 3 --- examples/fluids/src/setuplibceed.c | 3 --- 2 files changed, 6 deletions(-) diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index 458bb2fb0d..a52e3921b0 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -314,9 +314,6 @@ PetscErrorCode PrintRunInfo(User user, Physics phys_ctx, ProblemData *problem, M // ----------------------------------------------------------------------------- // libCEED functions // ----------------------------------------------------------------------------- -// Utility function - essential BC dofs are encoded in closure indices as -(i+1). -PetscInt Involute(PetscInt i); - // Utility function to create local CEED restriction PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt label_value, PetscInt dm_field, CeedElemRestriction *elem_restr); diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index ea028b05cb..c658f87a98 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -15,9 +15,6 @@ #include -// Utility function - essential BC dofs are encoded in closure indices as -(i+1). -PetscInt Involute(PetscInt i) { return i >= 0 ? i : -(i + 1); } - // Utility function to create local CEED restriction PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt label_value, PetscInt dm_field, CeedElemRestriction *elem_restr) { From 512ffa07b63531f745ea1b401cb5eaff13f61163 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 25 Jul 2023 11:22:25 -0600 Subject: [PATCH 04/15] fluids: Set quadrature points for PetscFE --- examples/fluids/navierstokes.c | 2 +- examples/fluids/navierstokes.h | 2 +- examples/fluids/src/setupdm.c | 22 +++++++++++----------- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/examples/fluids/navierstokes.c b/examples/fluids/navierstokes.c index 704f0b9b40..ec5eecbeb2 100644 --- a/examples/fluids/navierstokes.c +++ b/examples/fluids/navierstokes.c @@ -161,7 +161,7 @@ int main(int argc, char **argv) { } // -- Set up DM - PetscCall(SetUpDM(dm, problem, app_ctx->degree, bc, phys_ctx)); + PetscCall(SetUpDM(dm, problem, app_ctx->degree, app_ctx->q_extra, bc, phys_ctx)); // -- Refine DM for high-order viz if (app_ctx->viz_refine) PetscCall(VizRefineDM(dm, user, problem, bc, phys_ctx)); diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index a52e3921b0..3ac6df7613 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -358,7 +358,7 @@ PetscErrorCode UpdateBoundaryValues(User user, Vec Q_loc, PetscReal t); PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, MatType, VecType, DM *dm); // Set up DM -PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc, Physics phys); +PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_extra, SimpleBC bc, Physics phys); // Refine DM for high-order viz PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, SimpleBC bc, Physics phys); diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index ba200a3a4a..c0c1a543d2 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -39,23 +39,23 @@ PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, MatType mat_type, V } // Setup DM -PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc, Physics phys) { +PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_extra, SimpleBC bc, Physics phys) { + PetscInt q_order = degree + q_extra; PetscFunctionBeginUser; { PetscBool is_simplex = PETSC_TRUE; - - // Check if simplex or tensor-product mesh + + // Check if simplex or tensor-product mesh PetscCall(DMPlexIsSimplex(dm, &is_simplex)); // Configure the finite element space and boundary conditions PetscFE fe; PetscInt num_comp_q = 5; DMLabel label; - PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, is_simplex, degree, PETSC_DECIDE, &fe)); + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, is_simplex, degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Q")); PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); PetscCall(DMCreateDS(dm)); - // Project coordinates to enrich quadrature space - { + { // Project coordinates to enrich quadrature space DM dm_coord; PetscDS ds_coord; PetscFE fe_coord_current, fe_coord_new; @@ -64,14 +64,14 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc PetscCall(DMGetCoordinateDM(dm, &dm_coord)); PetscCall(DMGetCoordinateDim(dm, &num_comp_coord)); - PetscCall(DMGetRegionDS(dm_coord, NULL, NULL, &ds_coord,NULL)); + PetscCall(DMGetRegionDS(dm_coord, NULL, NULL, &ds_coord, NULL)); PetscCall(PetscDSGetDiscretization(ds_coord, 0, (PetscObject *)&fe_coord_current)); PetscCall(PetscFEGetDualSpace(fe_coord_current, &fe_coord_dual_space)); PetscCall(PetscDualSpaceGetOrder(fe_coord_dual_space, &fe_coord_order)); - // Create FE for coordinates - fe_coord_order = 1 ; - PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_coord, is_simplex, fe_coord_order, degree, &fe_coord_new)); + // Create FE for coordinates + fe_coord_order = 1; + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_coord, is_simplex, fe_coord_order, q_order, &fe_coord_new)); PetscCall(DMProjectCoordinates(dm, fe_coord_new)); PetscCall(PetscFEDestroy(&fe_coord_new)); } @@ -157,7 +157,7 @@ PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, SimpleBC bc, if (i + 1 == user->app_ctx->viz_refine) d = 1; PetscCall(DMGetVecType(dm, &vec_type)); PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type)); - PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, bc, phys)); + PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, 0, bc, phys)); PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL)); if (!i) user->interp_viz = interp_next; else { From 586d684e51f47364a98f570e184010edd5000e09 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 25 Jul 2023 12:21:27 -0600 Subject: [PATCH 05/15] fluids: Fix GetRestrictionForDomain, Use CreateBasisFromPlex --- examples/fluids/navierstokes.h | 5 +- examples/fluids/problems/sgs_dd_model.c | 12 ++--- examples/fluids/src/differential_filter.c | 9 ++-- examples/fluids/src/grid_anisotropy_tensor.c | 19 +++---- examples/fluids/src/setuplibceed.c | 53 ++++++++----------- .../fluids/src/strong_boundary_conditions.c | 4 +- examples/fluids/src/turb_spanstats.c | 17 +++--- .../fluids/src/velocity_gradient_projection.c | 9 ++-- 8 files changed, 53 insertions(+), 75 deletions(-) diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index 3ac6df7613..53f9c38667 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -319,10 +319,13 @@ PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLab CeedElemRestriction *elem_restr); // Utility function to get Ceed Restriction for each domain -PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, CeedInt Q_dim, +PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, CeedElemRestriction *elem_restr_qd_i); +PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, + CeedQuadMode quadMode, CeedBasis *basis); + // Utility function to create CEED Composite Operator for the entire domain PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, CeedData ceed_data, Physics phys, CeedOperator op_apply_vol, CeedOperator op_apply_ijacobian_vol, CeedInt height, CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur, diff --git a/examples/fluids/problems/sgs_dd_model.c b/examples/fluids/problems/sgs_dd_model.c index 1e76b3dad4..4fd80f1f9c 100644 --- a/examples/fluids/problems/sgs_dd_model.c +++ b/examples/fluids/problems/sgs_dd_model.c @@ -66,7 +66,7 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c SGS_DD_Data sgs_dd_data = user->sgs_dd_data; CeedQFunction qf_multiplicity, qf_sgs_dd_nodal; CeedOperator op_multiplicity, op_sgs_dd_nodal; - CeedInt num_elem, elem_size, num_comp_q, num_qpts_1d, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso; + CeedInt num_elem, elem_size, num_comp_q, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso; PetscInt dim; CeedVector multiplicity, inv_multiplicity; CeedElemRestriction elem_restr_inv_multiplicity, elem_restr_grad_velo, elem_restr_sgs; @@ -78,7 +78,6 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c CeedElemRestrictionGetNumComponents(sgs_dd_setup_data->elem_restr_grid_aniso, &num_comp_grid_aniso); CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &num_elem); CeedElemRestrictionGetElementSize(ceed_data->elem_restr_q, &elem_size); - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); { // Get velocity gradient information CeedOperatorField op_field; @@ -86,8 +85,7 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_grad_velo); CeedElemRestrictionGetNumComponents(elem_restr_grad_velo, &num_comp_grad_velo); } -// 0 added FIXME - PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, num_qpts_1d, 0,0, &elem_restr_sgs, NULL, NULL)); + PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, -1, 0, &elem_restr_sgs, NULL, NULL)); CeedElemRestrictionCreateVector(elem_restr_sgs, &sgs_dd_data->sgs_nodal_ceed, NULL); // -- Create inverse multiplicity for correcting nodal assembly @@ -160,7 +158,7 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c // @brief Create CeedOperator to compute SGS contribution to the residual PetscErrorCode SGS_ModelSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_data, SGS_DD_ModelSetupData sgs_dd_setup_data) { SGS_DD_Data sgs_dd_data = user->sgs_dd_data; - CeedInt num_comp_q, num_comp_qd, num_comp_x, num_qpts_1d, num_nodes_1d; + CeedInt num_comp_q, num_comp_qd, num_comp_x; PetscInt dim; CeedQFunction qf_sgs_apply; CeedOperator op_sgs_apply; @@ -171,10 +169,8 @@ PetscErrorCode SGS_ModelSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_ CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x); - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); - CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); - CeedBasisCreateTensorH1Lagrange(ceed, dim, sgs_dd_data->num_comp_sgs, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_sgs); + PetscCall(CreateBasisFromPlex(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, CEED_GAUSS, &basis_sgs)); switch (user->phys->state_var) { case STATEVAR_PRIMITIVE: diff --git a/examples/fluids/src/differential_filter.c b/examples/fluids/src/differential_filter.c index b0970c487c..baf4cfbfe2 100644 --- a/examples/fluids/src/differential_filter.c +++ b/examples/fluids/src/differential_filter.c @@ -17,7 +17,7 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData ceed_data, CeedQFunctionContext diff_filter_qfctx) { DiffFilterData diff_filter = user->diff_filter; DM dm_filter = diff_filter->dm_filter; - CeedInt num_comp_q, num_comp_qd, num_qpts_1d, num_nodes_1d, num_comp_x; + CeedInt num_comp_q, num_comp_qd, num_comp_x; PetscInt dim; PetscFunctionBeginUser; @@ -25,8 +25,6 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd); - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); - CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); { // -- Create RHS MatopApplyContext CeedQFunction qf_rhs; @@ -64,9 +62,8 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData char field_name[PETSC_MAX_PATH_LEN]; CeedElemRestriction elem_restr_filter; CeedBasis basis_filter; -// 0 added FIXME - PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, num_qpts_1d, 0,0, &elem_restr_filter, NULL, NULL)); - CeedBasisCreateTensorH1Lagrange(ceed, dim, diff_filter->num_field_components[i], num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_filter); + PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, -1, 0, &elem_restr_filter, NULL, NULL)); + PetscCall(CreateBasisFromPlex(ceed, dm_filter, 0, 0, 0, i, CEED_GAUSS, &basis_filter)); PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i)); CeedOperatorSetField(op_rhs, field_name, elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); diff --git a/examples/fluids/src/grid_anisotropy_tensor.c b/examples/fluids/src/grid_anisotropy_tensor.c index f8c250a4b5..d73d7e49e5 100644 --- a/examples/fluids/src/grid_anisotropy_tensor.c +++ b/examples/fluids/src/grid_anisotropy_tensor.c @@ -19,7 +19,7 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce CeedQFunction qf_rhs_assemble, qf_mass; CeedBasis basis_grid_aniso; PetscInt dim; - CeedInt num_qpts_1d, num_nodes_1d, q_data_size; + CeedInt q_data_size; MPI_Comm comm = PetscObjectComm((PetscObject)user->dm); KSP ksp; @@ -55,14 +55,10 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce } // -- Get Pre-requisite things - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); - CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); -//HACK to get to compile but won't work until code changed to get basis from PetscFE added extra 0 to arg list FIXME - PetscCall( - GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, num_qpts_1d, 0,grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); - CeedBasisCreateTensorH1Lagrange(ceed, dim, grid_aniso_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grid_aniso); + PetscCall(GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, -1, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); + PetscCall(CreateBasisFromPlex(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, CEED_GAUSS, &basis_grid_aniso)); // -- Build RHS operator CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorProjection, AnisotropyTensorProjection_loc, &qf_rhs_assemble); @@ -137,17 +133,15 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso, CeedVector *aniso_colloc_ceed, PetscInt *num_comp_aniso) { - CeedInt q_data_size, num_qpts_1d, num_qpts; + CeedInt q_data_size, num_nodes; CeedQFunction qf_colloc; CeedOperator op_colloc; PetscFunctionBeginUser; - // -- Get Pre-requisite things *num_comp_aniso = 7; - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); + CeedBasisGetNumNodes(ceed_data->basis_q, &num_nodes); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); -// here too FIXME - PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_qpts_1d, 0,*num_comp_aniso, NULL, NULL, elem_restr_grid_aniso)); + PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_nodes, *num_comp_aniso, NULL, NULL, elem_restr_grid_aniso)); // -- Build collocation operator CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc); @@ -157,7 +151,6 @@ PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User use CeedOperatorCreate(ceed, qf_colloc, NULL, NULL, &op_colloc); CeedOperatorSetField(op_colloc, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); CeedOperatorSetField(op_colloc, "v", *elem_restr_grid_aniso, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); - CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts); CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, aniso_colloc_ceed, NULL); diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index c658f87a98..fa50d10b7b 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -11,9 +11,8 @@ #include #include +#include #include "../navierstokes.h" -#include - // Utility function to create local CEED restriction PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt label_value, PetscInt dm_field, @@ -34,10 +33,9 @@ PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLab // Utility function to get Ceed Restriction for each domain PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, - CeedInt Q_dim, CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, + CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, CeedElemRestriction *elem_restr_qd_i) { DM dm_coord; -// now passed in/out CeedInt Q_dim CeedInt loc_num_elem; PetscInt dim; CeedElemRestriction elem_restr_tmp; @@ -45,7 +43,6 @@ PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel PetscCall(DMGetDimension(dm, &dim)); dim -= height; -// now passed in/out for simplex set by PetscFE Q_dim = CeedIntPow(Q, dim); PetscCall(CreateRestrictionFromPlex(ceed, dm, height, domain_label, label_value, dm_field, &elem_restr_tmp)); if (elem_restr_q) *elem_restr_q = elem_restr_tmp; if (elem_restr_x) { @@ -58,8 +55,7 @@ PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel } if (elem_restr_qd_i) { CeedElemRestrictionGetNumElements(elem_restr_tmp, &loc_num_elem); - CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q_dim, q_data_size, q_data_size * loc_num_elem * Q_dim, CEED_STRIDES_BACKEND, - elem_restr_qd_i); + CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q, q_data_size, q_data_size * loc_num_elem * Q, CEED_STRIDES_BACKEND, elem_restr_qd_i); } if (!elem_restr_q) CeedElemRestrictionDestroy(&elem_restr_tmp); PetscFunctionReturn(PETSC_SUCCESS); @@ -78,11 +74,12 @@ PetscErrorCode AddBCSubOperator(Ceed ceed, DM dm, CeedData ceed_data, DMLabel do CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); // ---- CEED Restriction - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, Q_sur, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, - &elem_restr_qd_i_sur)); + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, + &elem_restr_x_sur, &elem_restr_qd_i_sur)); if (jac_data_size_sur > 0) { // State-dependent data will be passed from residual to Jacobian. This will be collocated. - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, Q_sur, num_qpts_sur, jac_data_size_sur, NULL, NULL, &elem_restr_jd_i_sur)); + PetscCall( + GetRestrictionForDomain(ceed, dm, height, domain_label, label_value, 0, num_qpts_sur, jac_data_size_sur, NULL, NULL, &elem_restr_jd_i_sur)); CeedElemRestrictionCreateVector(elem_restr_jd_i_sur, &jac_data_sur, NULL); } else { elem_restr_jd_i_sur = NULL; @@ -213,7 +210,6 @@ PetscErrorCode SetupBCQFunctions(Ceed ceed, PetscInt dim_sur, PetscInt num_comp_ PetscFunctionReturn(PETSC_SUCCESS); } - // ----------------------------------------------------------------------------- // Convert DM field to DS field // ----------------------------------------------------------------------------- @@ -233,16 +229,15 @@ PetscErrorCode DMFieldToDSField(DM dm, DMLabel domain_label, PetscInt dm_field, if (dm_field == fields[i]) { *ds_field = i; break; - } + } } PetscCall(ISRestoreIndices(field_is, &fields)); if (*ds_field == -1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); } - // ----------------------------------------------------------------------------- // Utility function - convert from DMPolytopeType to CeedElemTopology // ----------------------------------------------------------------------------- @@ -261,7 +256,6 @@ CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { } } - // ----------------------------------------------------------------------------- // Create libCEED Basis from PetscTabulation // ----------------------------------------------------------------------------- @@ -314,8 +308,7 @@ PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, } } - // Convert to libCEED orientation - { + { // Convert to libCEED orientation PetscBool is_simplex = PETSC_FALSE; IS permutation = NULL; const PetscInt *permutation_indices; @@ -356,15 +349,14 @@ PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, PetscCall(PetscFree(interp)); PetscCall(PetscFree(grad)); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); } - // ----------------------------------------------------------------------------- // Get CEED Basis from DMPlex // ----------------------------------------------------------------------------- -PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedQuadMode quadMode, - CeedBasis *basis) { +PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, + CeedQuadMode quadMode, CeedBasis *basis) { PetscDS ds; PetscFE fe; PetscQuadrature quadrature; @@ -408,10 +400,9 @@ PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedI CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, quadMode, basis); } - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); } - PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { PetscFunctionBeginUser; @@ -419,32 +410,30 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App // Set up CEED objects for the interior domain (volume) // ***************************************************************************** const PetscInt num_comp_q = 5; - const CeedInt dim = problem->dim, num_comp_x = problem->dim, q_data_size_vol = problem->q_data_size_vol, jac_data_size_vol = num_comp_q + 6 + 3, - P = app_ctx->degree + 1, Q = P + app_ctx->q_extra; + const CeedInt dim = problem->dim, num_comp_x = problem->dim, q_data_size_vol = problem->q_data_size_vol, jac_data_size_vol = num_comp_q + 6 + 3; CeedElemRestriction elem_restr_jd_i; CeedVector jac_data; - CeedInt num_qpts; + CeedInt num_qpts; // ----------------------------------------------------------------------------- // CEED Bases // ----------------------------------------------------------------------------- - DM dm_coord; + DM dm_coord; PetscCall(DMGetCoordinateDM(dm, &dm_coord)); PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, CEED_GAUSS, &ceed_data->basis_q)); PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, CEED_GAUSS, &ceed_data->basis_x)); PetscCall(CeedBasisCreateProjection(ceed_data->basis_x, ceed_data->basis_q, &ceed_data->basis_xc)); - // --- Get number of quadrature points for the boundaries CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts); // ----------------------------------------------------------------------------- // CEED Restrictions // ----------------------------------------------------------------------------- // -- Create restriction - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, num_qpts, q_data_size_vol, &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, num_qpts, q_data_size_vol, &ceed_data->elem_restr_q, &ceed_data->elem_restr_x, &ceed_data->elem_restr_qd_i)); - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, num_qpts, jac_data_size_vol, NULL, NULL, &elem_restr_jd_i)); + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, num_qpts, jac_data_size_vol, NULL, NULL, &elem_restr_jd_i)); // -- Create E vectors CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_ceed, NULL); CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &user->q_dot_ceed, NULL); @@ -609,9 +598,9 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App // CEED Bases // ----------------------------------------------------------------------------- - DMLabel label = 0; + DMLabel label = 0; PetscInt face_id = 0; - PetscInt field = 0; // Still want the normal, default field + PetscInt field = 0; // Still want the normal, default field PetscCall(CreateBasisFromPlex(ceed, dm, label, face_id, height, field, CEED_GAUSS, &ceed_data->basis_q_sur)); PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, face_id, height, field, CEED_GAUSS, &ceed_data->basis_x_sur)); PetscCall(CeedBasisCreateProjection(ceed_data->basis_x_sur, ceed_data->basis_q_sur, &ceed_data->basis_xc_sur)); diff --git a/examples/fluids/src/strong_boundary_conditions.c b/examples/fluids/src/strong_boundary_conditions.c index 10aabc0e97..a2ee39d6c2 100644 --- a/examples/fluids/src/strong_boundary_conditions.c +++ b/examples/fluids/src/strong_boundary_conditions.c @@ -46,8 +46,8 @@ PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, CeedData ceed_data, DM dm, Problem // Compute contribution on each boundary face for (CeedInt i = 0; i < bc->num_inflow; i++) { // -- Restrictions - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, Q_sur, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, &elem_restr_x_sur, - &elem_restr_qd_sur)); + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, + &elem_restr_x_sur, &elem_restr_qd_sur)); CeedElemRestrictionCreateVector(elem_restr_q_sur, &multiplicity, NULL); CeedElemRestrictionGetMultiplicity(elem_restr_q_sur, multiplicity); CeedElemRestrictionGetNumElements(elem_restr_q_sur, &num_elem); diff --git a/examples/fluids/src/turb_spanstats.c b/examples/fluids/src/turb_spanstats.c index 86150c2086..798ffd4991 100644 --- a/examples/fluids/src/turb_spanstats.c +++ b/examples/fluids/src/turb_spanstats.c @@ -178,24 +178,27 @@ PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, CeedElemRestriction elem_re PetscErrorCode SpanStatsSetupDataCreate(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem, SpanStatsSetupData *stats_data) { DM dm = user->spanstats.dm; PetscInt dim; - CeedInt P, Q, num_comp_x, num_comp_stats = user->spanstats.num_comp_stats; + CeedInt num_qpts_child, num_comp_x, num_comp_stats = user->spanstats.num_comp_stats; Vec X_loc; PetscFunctionBeginUser; PetscCall(PetscNew(stats_data)); PetscCall(DMGetDimension(dm, &dim)); - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &Q); - CeedBasisGetNumNodes1D(ceed_data->basis_q, &P); -// inserted this zero just to get it to compile FIXME - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, Q, 0, problem->q_data_size_sur, &(*stats_data)->elem_restr_parent_stats, + CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_child); + num_qpts_child *= dim; + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, num_qpts_child, problem->q_data_size_sur, &(*stats_data)->elem_restr_parent_stats, &(*stats_data)->elem_restr_parent_x, &(*stats_data)->elem_restr_parent_qd)); CeedElemRestrictionGetNumComponents((*stats_data)->elem_restr_parent_x, &num_comp_x); CeedElemRestrictionCreateVector((*stats_data)->elem_restr_parent_x, &(*stats_data)->x_coord, NULL); CeedElemRestrictionCreateVector((*stats_data)->elem_restr_parent_qd, &(*stats_data)->q_data, NULL); - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, &(*stats_data)->basis_x); - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_stats, P, Q, CEED_GAUSS, &(*stats_data)->basis_stats); + { + DM dm_coord; + PetscCall(DMGetCoordinateDM(dm, &dm_coord)); + PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, CEED_GAUSS, &(*stats_data)->basis_x)); + PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, CEED_GAUSS, &(*stats_data)->basis_stats)); + } PetscCall(CreateElemRestrColloc(ceed, num_comp_stats, (*stats_data)->basis_stats, (*stats_data)->elem_restr_parent_stats, &(*stats_data)->elem_restr_parent_colloc)); diff --git a/examples/fluids/src/velocity_gradient_projection.c b/examples/fluids/src/velocity_gradient_projection.c index 2b038bc6b8..649ca838f2 100644 --- a/examples/fluids/src/velocity_gradient_projection.c +++ b/examples/fluids/src/velocity_gradient_projection.c @@ -54,7 +54,7 @@ PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, User user, CeedData ce CeedBasis basis_grad_velo; CeedElemRestriction elem_restr_grad_velo; PetscInt dim; - CeedInt num_comp_x, num_comp_q, q_data_size, num_qpts_1d, num_nodes_1d; + CeedInt num_comp_x, num_comp_q, q_data_size; PetscFunctionBeginUser; PetscCall(PetscNew(&user->grad_velo_proj)); @@ -66,13 +66,10 @@ PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, User user, CeedData ce PetscCall(DMGetDimension(grad_velo_proj->dm, &dim)); CeedBasisGetNumComponents(ceed_data->basis_x, &num_comp_x); CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q); - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); - CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); -// inserted zero just to compile FIXME - PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, num_qpts_1d, 0, q_data_size, &elem_restr_grad_velo, NULL, NULL)); + PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, -1, 0, &elem_restr_grad_velo, NULL, NULL)); - CeedBasisCreateTensorH1Lagrange(ceed, dim, grad_velo_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grad_velo); + PetscCall(CreateBasisFromPlex(ceed, grad_velo_proj->dm, 0, 0, 0, 0, CEED_GAUSS, &basis_grad_velo)); // -- Build RHS operator switch (user->phys->state_var) { From bb735d7e5a307c22966f1520c76554bfbf0a6df8 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 25 Jul 2023 12:25:50 -0600 Subject: [PATCH 06/15] fluids: Remove Q_sur arg from SetupStrongBC_Ceed --- examples/fluids/navierstokes.h | 3 +-- examples/fluids/src/setuplibceed.c | 2 +- examples/fluids/src/strong_boundary_conditions.c | 9 ++++----- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index 53f9c38667..e131c3332e 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -445,8 +445,7 @@ PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User use // ----------------------------------------------------------------------------- // Setup StrongBCs that use QFunctions -PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt Q_sur, - CeedInt q_data_size_sur); +PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt q_data_size_sur); PetscErrorCode FreestreamBCSetup(ProblemData *problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference); PetscErrorCode OutflowBCSetup(ProblemData *problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference); diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index fa50d10b7b..fff6f8335d 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -643,7 +643,7 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App if (user->op_ijacobian) { CeedOperatorGetContextFieldLabel(user->op_ijacobian, "ijacobian time shift", &user->phys->ijacobian_time_shift_label); } - if (problem->use_strong_bc_ceed) PetscCall(SetupStrongBC_Ceed(ceed, ceed_data, dm, user, problem, bc, Q_sur, q_data_size_sur)); + if (problem->use_strong_bc_ceed) PetscCall(SetupStrongBC_Ceed(ceed, ceed_data, dm, user, problem, bc, q_data_size_sur)); if (app_ctx->sgs_model_type == SGS_MODEL_DATA_DRIVEN) PetscCall(SGS_DD_ModelSetup(ceed, user, ceed_data, problem)); } diff --git a/examples/fluids/src/strong_boundary_conditions.c b/examples/fluids/src/strong_boundary_conditions.c index a2ee39d6c2..9ecf6571b0 100644 --- a/examples/fluids/src/strong_boundary_conditions.c +++ b/examples/fluids/src/strong_boundary_conditions.c @@ -13,8 +13,8 @@ #include "../navierstokes.h" #include "../problems/stg_shur14.h" -PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, CeedData ceed_data, DM dm, ProblemData *problem, SimpleBC bc, Physics phys, CeedInt Q_sur, - CeedInt q_data_size_sur, CeedOperator op_strong_bc) { +PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, CeedData ceed_data, DM dm, ProblemData *problem, SimpleBC bc, Physics phys, CeedInt q_data_size_sur, + CeedOperator op_strong_bc) { CeedInt num_comp_x = problem->dim, num_comp_q = 5, num_elem, elem_size, stg_data_size = 1; CeedVector multiplicity, x_stored, scale_stored, q_data_sur, stg_data; CeedBasis basis_x_to_q_sur; @@ -154,8 +154,7 @@ PetscErrorCode DMPlexInsertBoundaryValues_StrongBCCeed(DM dm, PetscBool insert_e PetscFunctionReturn(PETSC_SUCCESS); } -PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt Q_sur, - CeedInt q_data_size_sur) { +PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt q_data_size_sur) { CeedOperator op_strong_bc; PetscFunctionBeginUser; @@ -177,7 +176,7 @@ PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User use PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); if (use_strongstg) { - PetscCall(SetupStrongSTG_Ceed(ceed, ceed_data, dm, problem, bc, user->phys, Q_sur, q_data_size_sur, op_strong_bc)); + PetscCall(SetupStrongSTG_Ceed(ceed, ceed_data, dm, problem, bc, user->phys, q_data_size_sur, op_strong_bc)); } } From 3ab7f99218a190603ed4e5d048864ed985ebab85 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 25 Jul 2023 15:53:11 -0600 Subject: [PATCH 07/15] fluids: Remove extra SetClosurePermutationTensor --- examples/fluids/src/setupdm.c | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index c0c1a543d2..7decc94487 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -107,7 +107,6 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_ if (!is_simplex) { DM dm_coord; PetscCall(DMGetCoordinateDM(dm, &dm_coord)); - PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); } From 6614ffc51eee812e658c60c27972de63ac0eabff Mon Sep 17 00:00:00 2001 From: James Wright Date: Wed, 26 Jul 2023 14:23:00 -0600 Subject: [PATCH 08/15] fluids: Correct quadrature order to auxillary DMs - Should factor in q_extra --- examples/fluids/problems/sgs_dd_model.c | 10 ++++++---- examples/fluids/src/differential_filter.c | 11 ++++------- examples/fluids/src/grid_anisotropy_tensor.c | 3 ++- examples/fluids/src/velocity_gradient_projection.c | 3 ++- 4 files changed, 14 insertions(+), 13 deletions(-) diff --git a/examples/fluids/problems/sgs_dd_model.c b/examples/fluids/problems/sgs_dd_model.c index 4fd80f1f9c..df86c78412 100644 --- a/examples/fluids/problems/sgs_dd_model.c +++ b/examples/fluids/problems/sgs_dd_model.c @@ -29,19 +29,20 @@ PetscErrorCode SGS_DD_ModelSetupDataDestroy(SGS_DD_ModelSetupData sgs_dd_setup_d } // @brief Create DM for storing subgrid stress at nodes -PetscErrorCode SGS_DD_ModelCreateDM(DM dm_source, DM *dm_sgs, PetscInt degree, PetscInt *num_components) { +PetscErrorCode SGS_DD_ModelCreateDM(DM dm_source, DM *dm_sgs, PetscInt degree, PetscInt q_extra, PetscInt *num_components) { PetscFE fe; PetscSection section; PetscInt dim; PetscFunctionBeginUser; - *num_components = 6; + *num_components = 6; + PetscInt q_order = degree + q_extra; PetscCall(DMClone(dm_source, dm_sgs)); PetscCall(DMGetDimension(*dm_sgs, &dim)); PetscCall(PetscObjectSetName((PetscObject)*dm_sgs, "Subgrid Stress Projection")); - PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, *num_components, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, *num_components, PETSC_FALSE, degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Subgrid Stress Projection")); PetscCall(DMAddField(*dm_sgs, NULL, (PetscObject)fe)); PetscCall(DMCreateDS(*dm_sgs)); @@ -326,7 +327,8 @@ PetscErrorCode SGS_DD_ModelSetup(Ceed ceed, User user, CeedData ceed_data, Probl // -- Create DM for storing SGS tensor at nodes PetscCall(PetscNew(&user->sgs_dd_data)); - PetscCall(SGS_DD_ModelCreateDM(user->dm, &user->sgs_dd_data->dm_sgs, user->app_ctx->degree, &user->sgs_dd_data->num_comp_sgs)); + PetscCall( + SGS_DD_ModelCreateDM(user->dm, &user->sgs_dd_data->dm_sgs, user->app_ctx->degree, user->app_ctx->q_extra, &user->sgs_dd_data->num_comp_sgs)); PetscCall(PetscNew(&sgs_dd_setup_data)); diff --git a/examples/fluids/src/differential_filter.c b/examples/fluids/src/differential_filter.c index baf4cfbfe2..55c8734da0 100644 --- a/examples/fluids/src/differential_filter.c +++ b/examples/fluids/src/differential_filter.c @@ -184,7 +184,7 @@ PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data, { // Create DM for filtered quantities PetscFE fe; PetscSection section; - PetscInt dim; + PetscInt dim, q_order = user->app_ctx->degree + user->app_ctx->q_extra; PetscCall(DMClone(user->dm, &diff_filter->dm_filter)); PetscCall(DMGetDimension(diff_filter->dm_filter, &dim)); @@ -195,8 +195,7 @@ PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data, if (diff_filter->do_mms_test) { diff_filter->num_field_components[0] = 1; - PetscCall( - PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe)); + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering - MMS")); PetscCall(DMAddField(diff_filter->dm_filter, NULL, (PetscObject)fe)); PetscCall(PetscFEDestroy(&fe)); @@ -206,15 +205,13 @@ PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data, PetscCall(PetscSectionSetComponentName(section, 0, 0, "FilteredPhi")); } else { diff_filter->num_field_components[0] = DIFF_FILTER_STATE_NUM; - PetscCall( - PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe)); + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering - Primitive State Variables")); PetscCall(DMAddField(diff_filter->dm_filter, NULL, (PetscObject)fe)); PetscCall(PetscFEDestroy(&fe)); diff_filter->num_field_components[1] = DIFF_FILTER_VELOCITY_SQUARED_NUM; - PetscCall( - PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[1], PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe)); + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[1], PETSC_FALSE, user->app_ctx->degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering - Velocity Products")); PetscCall(DMAddField(diff_filter->dm_filter, NULL, (PetscObject)fe)); PetscCall(PetscFEDestroy(&fe)); diff --git a/examples/fluids/src/grid_anisotropy_tensor.c b/examples/fluids/src/grid_anisotropy_tensor.c index d73d7e49e5..90a027c19a 100644 --- a/examples/fluids/src/grid_anisotropy_tensor.c +++ b/examples/fluids/src/grid_anisotropy_tensor.c @@ -35,7 +35,8 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce { // -- Setup DM PetscFE fe; PetscSection section; - PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grid_aniso_proj->num_comp, PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe)); + PetscInt q_order = user->app_ctx->degree + user->app_ctx->q_extra; + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grid_aniso_proj->num_comp, PETSC_FALSE, user->app_ctx->degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Grid Anisotropy Tensor Projection")); PetscCall(DMAddField(grid_aniso_proj->dm, NULL, (PetscObject)fe)); PetscCall(DMCreateDS(grid_aniso_proj->dm)); diff --git a/examples/fluids/src/velocity_gradient_projection.c b/examples/fluids/src/velocity_gradient_projection.c index 649ca838f2..a363108d5d 100644 --- a/examples/fluids/src/velocity_gradient_projection.c +++ b/examples/fluids/src/velocity_gradient_projection.c @@ -25,7 +25,8 @@ PetscErrorCode VelocityGradientProjectionCreateDM(NodalProjectionData grad_velo_ PetscCall(DMGetDimension(grad_velo_proj->dm, &dim)); PetscCall(PetscObjectSetName((PetscObject)grad_velo_proj->dm, "Velocity Gradient Projection")); - PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grad_velo_proj->num_comp, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); + PetscInt q_order = user->app_ctx->degree + user->app_ctx->q_extra; + PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grad_velo_proj->num_comp, PETSC_FALSE, degree, q_order, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Velocity Gradient Projection")); PetscCall(DMAddField(grad_velo_proj->dm, NULL, (PetscObject)fe)); PetscCall(DMCreateDS(grad_velo_proj->dm)); From dcb23840c24f1eb6fa54cd2504c5ae691025ffcf Mon Sep 17 00:00:00 2001 From: James Wright Date: Wed, 26 Jul 2023 14:23:28 -0600 Subject: [PATCH 09/15] fluids: Fix parent Restriction creation for spanstats --- examples/fluids/src/turb_spanstats.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/fluids/src/turb_spanstats.c b/examples/fluids/src/turb_spanstats.c index 798ffd4991..07b5bef7fb 100644 --- a/examples/fluids/src/turb_spanstats.c +++ b/examples/fluids/src/turb_spanstats.c @@ -178,16 +178,16 @@ PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, CeedElemRestriction elem_re PetscErrorCode SpanStatsSetupDataCreate(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem, SpanStatsSetupData *stats_data) { DM dm = user->spanstats.dm; PetscInt dim; - CeedInt num_qpts_child, num_comp_x, num_comp_stats = user->spanstats.num_comp_stats; + CeedInt num_qpts_child1d, num_comp_x, num_comp_stats = user->spanstats.num_comp_stats; Vec X_loc; PetscFunctionBeginUser; PetscCall(PetscNew(stats_data)); PetscCall(DMGetDimension(dm, &dim)); - CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_child); - num_qpts_child *= dim; - PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, num_qpts_child, problem->q_data_size_sur, &(*stats_data)->elem_restr_parent_stats, + CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_child1d); + CeedInt num_qpts_parent = CeedIntPow(num_qpts_child1d, dim); + PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, 0, num_qpts_parent, problem->q_data_size_sur, &(*stats_data)->elem_restr_parent_stats, &(*stats_data)->elem_restr_parent_x, &(*stats_data)->elem_restr_parent_qd)); CeedElemRestrictionGetNumComponents((*stats_data)->elem_restr_parent_x, &num_comp_x); CeedElemRestrictionCreateVector((*stats_data)->elem_restr_parent_x, &(*stats_data)->x_coord, NULL); From d08e812cf523dbe270072f583be9b0257896477b Mon Sep 17 00:00:00 2001 From: James Wright Date: Fri, 28 Jul 2023 11:03:58 -0600 Subject: [PATCH 10/15] fluids: Correct qdata creation for strong bcs --- examples/fluids/src/strong_boundary_conditions.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/fluids/src/strong_boundary_conditions.c b/examples/fluids/src/strong_boundary_conditions.c index 9ecf6571b0..8cd11b281e 100644 --- a/examples/fluids/src/strong_boundary_conditions.c +++ b/examples/fluids/src/strong_boundary_conditions.c @@ -46,12 +46,13 @@ PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, CeedData ceed_data, DM dm, Problem // Compute contribution on each boundary face for (CeedInt i = 0; i < bc->num_inflow; i++) { // -- Restrictions - PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, num_qpts_sur, q_data_size_sur, &elem_restr_q_sur, - &elem_restr_x_sur, &elem_restr_qd_sur)); + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, -1, -1, &elem_restr_q_sur, &elem_restr_x_sur, NULL)); CeedElemRestrictionCreateVector(elem_restr_q_sur, &multiplicity, NULL); CeedElemRestrictionGetMultiplicity(elem_restr_q_sur, multiplicity); CeedElemRestrictionGetNumElements(elem_restr_q_sur, &num_elem); CeedElemRestrictionGetElementSize(elem_restr_q_sur, &elem_size); + PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, elem_size, q_data_size_sur, NULL, NULL, &elem_restr_qd_sur)); + CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, num_comp_x, num_elem * elem_size * num_comp_x, CEED_STRIDES_BACKEND, &elem_restr_x_stored); CeedElemRestrictionCreateVector(elem_restr_x_stored, &x_stored, NULL); From 203fa65bbcc44051646de2b458329ae8ce83283a Mon Sep 17 00:00:00 2001 From: James Wright Date: Fri, 28 Jul 2023 16:14:50 -0600 Subject: [PATCH 11/15] fluids: Set tensor closer only for tensor-product elements Co-authored-by: Jeremy L Thompson --- examples/fluids/src/setupdm.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index 7decc94487..93e27c37ac 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -103,10 +103,10 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_ if (use_strongstg) PetscCall(SetupStrongSTG(dm, bc, problem, phys)); } - PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); if (!is_simplex) { DM dm_coord; PetscCall(DMGetCoordinateDM(dm, &dm_coord)); + PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); } From 3b1543e51a779733de557f1d633e6d3847de03bc Mon Sep 17 00:00:00 2001 From: James Wright Date: Mon, 31 Jul 2023 11:41:02 -0600 Subject: [PATCH 12/15] fluids: Make shocktube test less brittle - Small, insignificant changes to the coordinate locations make a dramatic difference to the initial condition - For the test case, there is a node right at the midpoint of the domain. Adding the epsilon fudge factor makes the if statement more stable --- examples/fluids/qfunctions/shocktube.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluids/qfunctions/shocktube.h b/examples/fluids/qfunctions/shocktube.h index 1552cc9e6b..ff26144f79 100644 --- a/examples/fluids/qfunctions/shocktube.h +++ b/examples/fluids/qfunctions/shocktube.h @@ -94,7 +94,7 @@ CEED_QFUNCTION_HELPER CeedInt Exact_ShockTube(CeedInt dim, CeedScalar time, cons CeedScalar rho, P, u[3] = {0.}; // Initial Conditions - if (x <= mid_point) { + if (x <= mid_point + 200 * CEED_EPSILON) { rho = rho_high; P = P_high; } else { From 168d77717a8ff02ebc93e9c7175bd1c6afd3f594 Mon Sep 17 00:00:00 2001 From: James Wright Date: Mon, 31 Jul 2023 17:12:38 -0600 Subject: [PATCH 13/15] fluids: Fix SGS test - Change the test so that it passes on all hardware. It's a very sentitive test and the changes to the basis creation changed the initial conditions *just* enough to make the test not happy. Should really try and make this more robust at some point... --- examples/fluids/navierstokes.c | 2 +- ...s-navierstokes-blasius-sgs-data-driven.bin | Bin 5272 -> 5272 bytes 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluids/navierstokes.c b/examples/fluids/navierstokes.c index ec5eecbeb2..157485dc77 100644 --- a/examples/fluids/navierstokes.c +++ b/examples/fluids/navierstokes.c @@ -18,9 +18,9 @@ // ./navierstokes -ceed /cpu/self -problem density_current -degree 1 // ./navierstokes -ceed /gpu/cuda -problem advection -degree 1 // +//TESTARGS(name="blasius_SGS_DataDriven") -ceed {ceed_resource} -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 2e-9 -state_var primitive -ksp_rtol 1e-12 -snes_rtol 1e-12 -stg_mean_only -stg_fluctuating_IC -test_type solver -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin //TESTARGS(name="blasius_aniso_diff_filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 5e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_aniso_vandriest.bin -diff_filter_monitor -ts_max_steps 0 -state_var primitive -diff_filter_friction_length 1e-5 -diff_filter_wall_damping_function van_driest -diff_filter_ksp_rtol 1e-8 -diff_filter_grid_based_width -diff_filter_width_scaling 1,0.7,1 //TESTARGS(name="blasius_iso_diff_filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 2e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_iso.bin -diff_filter_monitor -ts_max_steps 0 -diff_filter_width_scaling 4.2e-5,4.2e-5,4.2e-5 -diff_filter_ksp_atol 1e-14 -diff_filter_ksp_rtol 1e-16 -//TESTARGS(name="blasius_SGS_DataDriven") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 1e-9 -compare_final_state_atol 2e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin -state_var primitive //TESTARGS(name="gaussianwave_idl") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 2e-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-IDL.bin -dm_plex_box_faces 5,5,1 -ts_max_steps 5 -idl_decay_time 2e-3 -idl_length 0.25 -idl_start 0 //TESTARGS(name="turb_spanstats") -ceed {ceed_resource} -test_type turb_spanstats -options_file examples/fluids/tests-output/stats_test.yaml -compare_final_state_atol 1E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-turb-spanstats-stats.bin //TESTARGS(name="blasius") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 2E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius.bin diff --git a/examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin b/examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin index de45762d35a915874b4489ad4d3d6372b47d44b0..5aee7101debbd93f9e45c9667dd377f1694accd0 100644 GIT binary patch literal 5272 zcmWldcOX@79LAB=uy<%Fl)YuWFZX6-&$dEk6&b%I3ZX@!5G5f*GAdc&$Y|MaNJtqK zD#X3+mH2({-_LozpZC1u`#$Gv_&uV|!otE~v?nxSR5f0YMKH{Ug@r8<(^Jazo!Sj? zz_Q7x_>eZH+j~p6eICRCzdsl6Y?{UlA(~O{+hiOly!entV|JJQ)~-F*(EOs;+BHA_!1#m0oMTlnJ?~bpwYUrVOE0e% zqBaq9v(OxiyCeFl`*IN*V6)Vy@hS&sVIOYp!~#pp ze?OklKt0WCZ9@Hvz)#ybo!kQTWz$abiC`h0P{ow|C`6HT%-pS<1G+?pyTnh>_0!I} z8)IqV)xABNUm&U=mxE79GU(qHRNtAv`snS>ZTX8OT%CSc%5)HnG=nydLuJ6Xy|PRb z0KP)8jEXpxh>l-U(L9HGI!fnaPCo(vLT$XA4a{$?D*Id#-Xx|UazcF*qPTiD8uOar zfR*=CU4k(Gk9=}_2U4+2jp|LKY!^h`lqr9!?-BS{N9%oI|34PnX-!z*t(8ThU%tr_ zjH}yZ4@*At?qOR-y-9Rajp-oZlPlMGzX$%jpIu}v*3LHP zVO^?6RNn(3Q9WMZ_jR1AY7D9lQF0k1y8j zc#kz_)CSe~KB5l~pZ6DpHh|7qr6PO;(;pV^=iA?a^-}Hy3t6Y4{zIoXWzBwt{+Jls zvcLy^jfZX?U$8;aYp1R0g@{_+m|>Yw1O9aPz=7SsAK+UEIEVF>cG?Ivwi9%<%v1Yx zt3VeGuDx^&{HFA~tNXEDO@$C^MGWd&H=-CaxCL~cYibL+n4a^dEBci=#yetq?L25` zz>a~Wi1j$&@3{RAv3&5KpBPyw!N#g8pHlmz5KY)%v!Npk_}_EMT|WW*SGR;)kBqaYKuaaJHf7@4a9Q54JExa!}d)ejRZ zHs@keSK(8E@En33+bn&2YKl4k4SOc*fN!`pwxAA^XmQsI!fes!`+YYH%et9-qhC?5 zzFE9qk{{b+lb15e!3kNYC5TnxvLf^R*{&u_zJUI!nKsdJ6_f9|66R8C5M8d)NAuNd z@c#~b^WSCQE3}#Tg=5l~^xlAZBqpB`B*W}~hbKGL>M>bt?Uqab1X2G- zoj>{8XTUF+*H9)5{JN92i-Fi=SjN>g^E{#%?C^cG-WB-P{)JqBp?+^8spfqyt)#{U_0;s|E& zXLJXZe}epI%2}`)0r~M=^{3md@ch1AmjYpbJ#nQDaVR&i-q47VeRm91-K9~+0{fwU zn)PI9wMYi=^?M_oVD+a?dg{lGRlvdlY)>9 zT-R14%L7f)zA5CFHV~_?6FX|oWa2>9S!!_z%zu9Fy`!VOXjU{hExV2t^-bQIkS`N~ z`j12W(;$Bo8t6P9ttDzqJ#R?ZeMCK)68+uxe}n$hIsI1>rWZ%EKK2#FYept`0@ZZU z!p0_|xph6%=Z+JSiGu#Cpk{Gf7yVtIvBfa00X^q#Qty{s8ls*z(D=BKnI9An%#BonuJK?eM-S#~EN`hAlOdKg9~jC{%Y%Qc zH9_eX=yCH(51wMKF^07IVmhL%-)=5ny9|ELf_k?rpzlgdbqK{ABd1qw&3;3^i0yQ# z(gy$eve#cu!TB*y|C7Usg$IT&`Rb>GOS;=UK)4-a*amW28B^opG8!U9L#&3PLj zKmHZ(ozhbS|J57i8N*QjlM^FA6bm%W4@}%HMLo+NIN_}|_{BDt-Ef0^@Ve)U?JO43 zP*>d$WsN9)x|T{U{NUG1tgnI4MM9ml~oTUTFM%Y>3I1udYl0em5;l z>MWL!UDe2Zo@)JgMK>jwG5Jqt+3342pefnenh=Ee#DWj z2mC2w#u|HoUn<5m`4exwFQsHnl^_@^VXd7mOTbUy+tLg9q@e4T&6!HP^~$fry?i?m z<>-DT<(FaLmpDi}7YhDP!?BS>yrZ+^TgT7_^uhkku$PrP@MG)RBq85Db`frX{kF?m zL_B_215vY4VP!M(`=yL^8`}i<^YaAz&tiFr$&e7zD#4iF;BfT-3HX;?saAr(r*u_D z9K~{!8%57ALw;;N{Cp_S4g6BHqZV*K=S%w>eP4u?qj(hV8kL}4Rn~Pplx5(v+d8<+ z06!+ar?CsGxxXet#7-jW@T)TI?vvoBtbg%01N&iWSc|07 z(kZr;VAQu?e4s=IzF*{7o6W%g!d`P}32To3ym3O79ra#qk7h8K`BpB+@YPz#uTHL& zk*!$oX2%n={nydQYa6ffZ4QP0^)yz!RT=mz+R3R0vB8zZ=Z2U5AX@bo@w(@bA6Qi! z>pWn8=KI$7c*kM=?ZW5EY?2AaSxb3guE)SX`d`AJ1n?&xH%>tPXKn65HWjFUvom>} z;#J@)CwJyUKYr-HM@;x5#@nAS@D7sDpr^N{FIxie*-f0ML*V;RH_OLz85=9a{~QbD zLW{>Y%uc#a0{_U^Ri_2uyX&?e62droUzq;Y9|U8k_5H}kI^a9jJm-ddRWP3|X|w}l z2Nzxr2~#xS@7*B%@C5L+3=LDxL;sp9aQ-KVO++nM?Nc()@IfC7BVKE`zfwNN=Y9vi z@7RL&AxzRQso>C_Lv(4gFI7z=z>n1N>#G62sh;zLa%}R#?QZ-w4uU?tQDt{fIq=Pj zscg)A#hvSV{wOxd)16Fq9Yllq-hnNm{J__LuW51*`1k2mkHWCY>*mkiF`Lj(+acjr z#a`f_o_O^SzF&(@$B+-!Ve+~;@ta&)XmOS@H7jr(_?NA%FR=jMXYh`H2qxuzd!sdM zK`j1F%R3Rn{JyhcaWMzIkbT|K<}6IAS{9YQ=#7TLPKaz^Q)KeJqE12n$eBwLKldA( zOcF+=T?f$*c>~%uzeV8tB`M9q|BE8W{({k7Oy2S!N7f?{Eh)r_n#+&1M$AZq*>RC8A2D)*J?Y-dmBi= zH}v}6d%RowQ<%hd@7X@PSoA&7P5tFYMc|)tRPru`{kN<=_!I7b(a$rj51J6o zb#LWLL?!U!n*;qePty7! zzdIWly@&6&W3h2>%%_?B(_+V(fuG-_dZZ7NrplgkSd9^jlS69>GJJoonwv>kvjKlw zPSbwKZ>t>Jht~W?L*=)Mq^>J4`AvHc&%^n%tL5T*i%C`w@~&KdjK;ZgJR^f`nDhT+ zjiD0o#q|QVYyv(lH{AvPU&xmHN)KKK{2Mc2Cn104FZd{)gnYJ^`j$`f5JBS;syw$K z%A9{PS1|+j|4Nkap8`x?FBhPxU5Y+-;2V~8%>HlG=bs}5{P7{z%Z8YA-bEodoEuGh z@NcN!6b^jQ)#xj|;WuvqirHvE4mprCdP`F~?Q>o`jj8YX@pk&zSxK04bzv>y1v^{=05 zW73{kKShyG#LQNw&ldvU0zatr&$%)G~nxA%{?r^ Ku zp4QJPSJ|yb4A7aR;u(3CCy+1Cd2${j@`q=fP0_5?qjDX_;#T23!?Jp(J`nSER{h-# zd3b0^fB_nnAGC07{-H+Am5>ahn`e^dg6c%yzs zby*IH8Gh2anBE2b&8*cUQN;Wv%Z3I_x9gu$2;ViN)c>COTl@s{*O~a79fC_+53zqLu+}E{g)U7mAXn2Wo!W{SZN5I~OY2A0(|l9sO)o?L64#QPHy@G{ndcL ztjmyxB;GOlhFK&2OIctmN7`3wBrC+9L7)D#Xm35_yNe3@He)W#M{SAsI_%U}yyGfxt2eRJy zaht}O3jLMv-W*@ZL!A;A4pF>QU;Nu7 z*DmBk3mSWsdoljFqv z9@7O&&kmx?jV_~2{prwOoS&Q;gFGzO_`qY#FEE<2ms10|O#3)r=Td?`A1 zr#u_MBEd->-P+d#H{|l|ujC%@>NfSnPWLagx(>6m;k{NKmU~VVVHx3;H(s~qP{xi=ZVicP&}aF}wHS(*N~;gJ2+(8sp_?D(4?aNUQcX41 zH-ms*Whr+GL7yg(PtnH;iHA7Fyv9)p=g$b&l@#ckO?SRXfZUYhlAR`&&xm@@Wg?A= zofw4csuO`f?fWWh3*>PLtgcV7{4V=}Bt-{Q#9cALtIz~}wwd(sEZ{3xbBgp~l{ZB) zavYQBL44_5wbc#a8=A54Q1Ma)PsilLBUm+LI3X>_2-O^Jxy2E+ANrp5eNs9wUwOj@a>OHwdr+_~$BQiM)d18K2*bF}Qc`9)A_fd2|$bGrY%NzQf7J`m@ z@KQ;KOx;E~#$6oi)uxrGxhwJ+OR53(SX);dwuALwAt$C+KZVsLrVA8?%@O6lQ*4xo z>(GxAAI?#R`I=av4#ycb8~c3Icy9@oKt1Y z$Vn85Jm_8e4R&GF847>Y5GXd zz(>?p=|Xp7%YEo)FRSIV!+eeA7jv;#W5PCRZeu5E%XP0U|6KNq`z+W3k$KCvBz{*D{2mau+GY&M@6xQa=RdG|bnjv{Jn<}vEL2=MDB?#LYk zz7g;C#Yjxv%pZHYD-KcDswT+t?9k8SWN|Yj_)b<@eVFvpcdl(}7Evi5mtOi1`^)R` z9XvlE&xtOJ^TMRQkKJxeyAieQhp!Lq5%hVRP}Lyt&tFh>4Z)&PzYcXk|``--%WzAnM7NJ%1wa68tss?XbQ|9A4dPv&i$r?yyBEM)A@hEZ{i*MYlqqde%Y*Pco2N!<=&Gp|Kyg4H+M1V zZvt=X2O4^MBb?v-S1iG2+h*SjdA62_vI8bF1;@52HKBIzf`iv;+y2j2ILroow@f9= zU`z^*wkr_+ih7i{A8)qa5Bvw2tdi!yKfxTjSdYmf{ufWbe}-C=ZMO!I{0M&f{qZd7AWQ*7tD!1s65qv{2jWZ@)skiTlsT z>F;t6`2Kz~o_jDUH13><4GVgHrRL6xnIiCO5~f&T{a2{mER89{WTAN8+*WVY=oEdV zQgQ?O`SP>Q`(VE5!saS1Oj?%aPxSOeJwn3u4~1_N{4NHCTad?pn7BNK$qdJr&3fOV zW)o^~cZ(pwm)4(y{i`H(M4^cdlOjjwsc|LfC5xG@MT->h?>!Cwkqvw?=0Vm=m@IJ1 z(Y=8M)&Fwzca{o=e!>2>d|h`2@ShTIrCK(GNg*D*fn(2G+eBYo)f=}W_)6nBuzyvk z`Rs_%#$@3GoS$di(A&L&&ZDxm1Yh*o<2;yeEPTl8CnhartP96*Q<>!TI-@X{WlNenuZ|y*0w$_=F`zydV!Myt2dz#|%@h4T^9a8f16I3;G)x z!ow}ZG5b&c4#?@WK9|0L6Gd_NXusT zh8QO?|L2(My++8x&~tk$@JEG?l)<2F{(v>5zTyys8S=h+C2NW)KKYp;`fQU-F^tUNmTdw^Wcvm zmo8=>L+3K5mACR-hyE0o=_g|SM^T;~xB&jxt<~tzf&>-0PYtshK!4NAhB$b?mj}0V ze9^;n{2a9r>J~^r;L?=yTs`#Hd28JXe+V9xi46sROrres&Op-9!SmeFbp-$B6Ln9> z)y{D249AS;nJTk=%#o#F28YM~aOl4pqLB&zO*qH4x*7a&q91>hMWmZWOe&OV=(Dw1 z&>x5VKg;KlHQoZ`s;a3;sCP$;dv4toGc^jm00Jzmo2)tOR+e$GH=`z#kQx_R8-> zW(|T5j|ap-pFzwxc$)A>{mBUh@W)g6s%t#RJMkZ9TU{UYv&-gq{_l^;$A^uWJUZZ?Q*$ct81gkU6eneoyYKeHJ3fO`_oy>v!;Q>8Gj&mzwou^H~8aX{RSlm1@`vieGaP7UwD{8 zCk1)L)*g?XW@qjwc+Dtn3;gff?q(A2qZl_6raJJ)?*XUZ7ozmD+J<^$5$LnbH}8b^dujEN z@gZDQNwUyf&tA*X7k>xpW2YpT>FSebqe{}4S^;idg6z5&eNkYxhZ(KW~^I`v3Pg62D3w$}ja+?6~ z#|?|Cm(~zPeP))c(F*!;tUN+HV7_L7QvetEqn+Kr3_rWAJLB{e%%IQx!BCqF`~Z?n&nEE4aMsz}x2S#0d|+~T0Qj9klwvO6 z>(mxSuYf|!iBo+o)pVJb|d&Mip*z#ub-jfXbAkI9NqTM=<#90FVhcg(bM0tb<#t`{eMl9)Ft-Mte;kk{NRtSVtzi6sQoaNygh<=e)-s1-@^An ziHG-@p;Yjf1BYpa`_OZRLq}}=n}A=-VsiE)@JSRd6AJjttT@wH9(p9lwWKF`1^W4$ zT5Pgl{@FmWA${=2KNk7<^r(xbn48`T-+vpg`#B2TAy4tFx~vZ8Z?rUdED=2(kCS;f zAxrSZDv#2DA8LLp6Z}1_yd=fc4!ziHDfm$A1n}=)8*PgNzP#i2E$f)PC!QyTDGN3D z^cQx-#}NEB%5@UVx5`w0xB&k6qeJ?&2Hz# zqNbL0?n{^X2tG+n3f6Z?I@M*@IZTR@P4sSmik{QSiVYisf3Dpv_+xhi_#%5%#;3p^ zPw4Usoj^^F>^JOg64$>VXv9eq=3l-b8ul9eu{gv1@MY9>GLKhO{3pTpJ^AMU`{Rr^ zUljP`x!C%lcc^Wfjb47L6TvsrR+R+)m8C_=7vS&IKl%zBtxW@a%42Q{1HY=Pqu&|$ z8t%ETdg1!NJbLwWG8#H9vZH`JMVvonx`4R9276k}v%w#C^DBtfBitWMafAO00>5iE jN(Ay`-Z?`Cxc;)!`?Fo)Y5!FMK>{C+GeL!~5Wj From 7d7afb295391ff73fa46eafb343a4e665119b143 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 1 Aug 2023 14:10:28 -0600 Subject: [PATCH 14/15] fluids: Address misc PR comments --- examples/fluids/navierstokes.h | 2 +- examples/fluids/src/setupdm.c | 11 +++++++---- examples/fluids/src/setuplibceed.c | 4 ++-- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index e131c3332e..edf042e66f 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -324,7 +324,7 @@ PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel CeedElemRestriction *elem_restr_qd_i); PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, - CeedQuadMode quadMode, CeedBasis *basis); + CeedQuadMode quad_mode, CeedBasis *basis); // Utility function to create CEED Composite Operator for the entire domain PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, CeedData ceed_data, Physics phys, CeedOperator op_apply_vol, diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index 93e27c37ac..8713f9eeb0 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -58,7 +58,7 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_ { // Project coordinates to enrich quadrature space DM dm_coord; PetscDS ds_coord; - PetscFE fe_coord_current, fe_coord_new; + PetscFE fe_coord_current, fe_coord_new, fe_coord_face_new; PetscDualSpace fe_coord_dual_space; PetscInt fe_coord_order, num_comp_coord; @@ -70,8 +70,10 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_ PetscCall(PetscDualSpaceGetOrder(fe_coord_dual_space, &fe_coord_order)); // Create FE for coordinates - fe_coord_order = 1; + PetscCheck(fe_coord_order == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER_INPUT, + "Only linear mesh geometry supported. Recieved %d order", fe_coord_order); PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_coord, is_simplex, fe_coord_order, q_order, &fe_coord_new)); + PetscCall(PetscFEGetHeightSubspace(fe_coord_new, 1, &fe_coord_face_new)); PetscCall(DMProjectCoordinates(dm, fe_coord_new)); PetscCall(PetscFEDestroy(&fe_coord_new)); } @@ -152,11 +154,12 @@ PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, SimpleBC bc, PetscCall(DMClearDS(dm_hierarchy[i + 1])); PetscCall(DMClearFields(dm_hierarchy[i + 1])); PetscCall(DMSetCoarseDM(dm_hierarchy[i + 1], dm_hierarchy[i])); - d = (d + 1) / 2; + d = (d + 1) / 2; + PetscInt q_order = d + user->app_ctx->q_extra; if (i + 1 == user->app_ctx->viz_refine) d = 1; PetscCall(DMGetVecType(dm, &vec_type)); PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type)); - PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, 0, bc, phys)); + PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, q_order, bc, phys)); PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL)); if (!i) user->interp_viz = interp_next; else { diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index fff6f8335d..dc8f3658b4 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -356,7 +356,7 @@ PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, // Get CEED Basis from DMPlex // ----------------------------------------------------------------------------- PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, - CeedQuadMode quadMode, CeedBasis *basis) { + CeedQuadMode quad_mode, CeedBasis *basis) { PetscDS ds; PetscFE fe; PetscQuadrature quadrature; @@ -397,7 +397,7 @@ PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedI CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim)); CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim)); - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, quadMode, basis); + CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, quad_mode, basis); } PetscFunctionReturn(PETSC_SUCCESS); From bf9fb9807cc3b51d328fd5cffe9d4d0f5039a0f5 Mon Sep 17 00:00:00 2001 From: James Wright Date: Tue, 1 Aug 2023 16:33:37 -0600 Subject: [PATCH 15/15] fluids: Remove quad_mode arg from CreateBasisFromPlex --- examples/fluids/navierstokes.h | 3 +-- examples/fluids/problems/sgs_dd_model.c | 2 +- examples/fluids/src/differential_filter.c | 2 +- examples/fluids/src/grid_anisotropy_tensor.c | 2 +- examples/fluids/src/setuplibceed.c | 13 ++++++------- examples/fluids/src/turb_spanstats.c | 4 ++-- examples/fluids/src/velocity_gradient_projection.c | 2 +- 7 files changed, 13 insertions(+), 15 deletions(-) diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index edf042e66f..b8086af5c9 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -323,8 +323,7 @@ PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, CeedElemRestriction *elem_restr_qd_i); -PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, - CeedQuadMode quad_mode, CeedBasis *basis); +PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedBasis *basis); // Utility function to create CEED Composite Operator for the entire domain PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, CeedData ceed_data, Physics phys, CeedOperator op_apply_vol, diff --git a/examples/fluids/problems/sgs_dd_model.c b/examples/fluids/problems/sgs_dd_model.c index df86c78412..99138a98fc 100644 --- a/examples/fluids/problems/sgs_dd_model.c +++ b/examples/fluids/problems/sgs_dd_model.c @@ -171,7 +171,7 @@ PetscErrorCode SGS_ModelSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_ CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd); CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x); - PetscCall(CreateBasisFromPlex(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, CEED_GAUSS, &basis_sgs)); + PetscCall(CreateBasisFromPlex(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, &basis_sgs)); switch (user->phys->state_var) { case STATEVAR_PRIMITIVE: diff --git a/examples/fluids/src/differential_filter.c b/examples/fluids/src/differential_filter.c index 55c8734da0..114dbc72b5 100644 --- a/examples/fluids/src/differential_filter.c +++ b/examples/fluids/src/differential_filter.c @@ -63,7 +63,7 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData CeedElemRestriction elem_restr_filter; CeedBasis basis_filter; PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, -1, 0, &elem_restr_filter, NULL, NULL)); - PetscCall(CreateBasisFromPlex(ceed, dm_filter, 0, 0, 0, i, CEED_GAUSS, &basis_filter)); + PetscCall(CreateBasisFromPlex(ceed, dm_filter, 0, 0, 0, i, &basis_filter)); PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i)); CeedOperatorSetField(op_rhs, field_name, elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); diff --git a/examples/fluids/src/grid_anisotropy_tensor.c b/examples/fluids/src/grid_anisotropy_tensor.c index 90a027c19a..2e19354629 100644 --- a/examples/fluids/src/grid_anisotropy_tensor.c +++ b/examples/fluids/src/grid_anisotropy_tensor.c @@ -59,7 +59,7 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); PetscCall(GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, -1, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); - PetscCall(CreateBasisFromPlex(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, CEED_GAUSS, &basis_grid_aniso)); + PetscCall(CreateBasisFromPlex(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, &basis_grid_aniso)); // -- Build RHS operator CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorProjection, AnisotropyTensorProjection_loc, &qf_rhs_assemble); diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index dc8f3658b4..4f763a1488 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -355,8 +355,7 @@ PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, // ----------------------------------------------------------------------------- // Get CEED Basis from DMPlex // ----------------------------------------------------------------------------- -PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, - CeedQuadMode quad_mode, CeedBasis *basis) { +PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedBasis *basis) { PetscDS ds; PetscFE fe; PetscQuadrature quadrature; @@ -397,7 +396,7 @@ PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedI CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim)); CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim)); - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, quad_mode, basis); + CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, CEED_GAUSS, basis); } PetscFunctionReturn(PETSC_SUCCESS); @@ -421,8 +420,8 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App DM dm_coord; PetscCall(DMGetCoordinateDM(dm, &dm_coord)); - PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, CEED_GAUSS, &ceed_data->basis_q)); - PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, CEED_GAUSS, &ceed_data->basis_x)); + PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, &ceed_data->basis_q)); + PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, &ceed_data->basis_x)); PetscCall(CeedBasisCreateProjection(ceed_data->basis_x, ceed_data->basis_q, &ceed_data->basis_xc)); CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts); @@ -601,8 +600,8 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App DMLabel label = 0; PetscInt face_id = 0; PetscInt field = 0; // Still want the normal, default field - PetscCall(CreateBasisFromPlex(ceed, dm, label, face_id, height, field, CEED_GAUSS, &ceed_data->basis_q_sur)); - PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, face_id, height, field, CEED_GAUSS, &ceed_data->basis_x_sur)); + PetscCall(CreateBasisFromPlex(ceed, dm, label, face_id, height, field, &ceed_data->basis_q_sur)); + PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, face_id, height, field, &ceed_data->basis_x_sur)); PetscCall(CeedBasisCreateProjection(ceed_data->basis_x_sur, ceed_data->basis_q_sur, &ceed_data->basis_xc_sur)); // ----------------------------------------------------------------------------- diff --git a/examples/fluids/src/turb_spanstats.c b/examples/fluids/src/turb_spanstats.c index 07b5bef7fb..634226c615 100644 --- a/examples/fluids/src/turb_spanstats.c +++ b/examples/fluids/src/turb_spanstats.c @@ -196,8 +196,8 @@ PetscErrorCode SpanStatsSetupDataCreate(Ceed ceed, User user, CeedData ceed_data { DM dm_coord; PetscCall(DMGetCoordinateDM(dm, &dm_coord)); - PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, CEED_GAUSS, &(*stats_data)->basis_x)); - PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, CEED_GAUSS, &(*stats_data)->basis_stats)); + PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, &(*stats_data)->basis_x)); + PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, &(*stats_data)->basis_stats)); } PetscCall(CreateElemRestrColloc(ceed, num_comp_stats, (*stats_data)->basis_stats, (*stats_data)->elem_restr_parent_stats, diff --git a/examples/fluids/src/velocity_gradient_projection.c b/examples/fluids/src/velocity_gradient_projection.c index a363108d5d..75730096e1 100644 --- a/examples/fluids/src/velocity_gradient_projection.c +++ b/examples/fluids/src/velocity_gradient_projection.c @@ -70,7 +70,7 @@ PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, User user, CeedData ce CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, -1, 0, &elem_restr_grad_velo, NULL, NULL)); - PetscCall(CreateBasisFromPlex(ceed, grad_velo_proj->dm, 0, 0, 0, 0, CEED_GAUSS, &basis_grad_velo)); + PetscCall(CreateBasisFromPlex(ceed, grad_velo_proj->dm, 0, 0, 0, 0, &basis_grad_velo)); // -- Build RHS operator switch (user->phys->state_var) {