diff --git a/examples/fluids/navierstokes.c b/examples/fluids/navierstokes.c index 704f0b9b40..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 @@ -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 30d96eeeed..b8086af5c9 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); @@ -326,6 +323,8 @@ 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, 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, @@ -361,7 +360,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); @@ -445,8 +444,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/problems/sgs_dd_model.c b/examples/fluids/problems/sgs_dd_model.c index 33eac19293..99138a98fc 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)); @@ -66,7 +67,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 +79,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 +86,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); } - - PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, num_qpts_1d, 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 +159,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 +170,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, &basis_sgs)); switch (user->phys->state_var) { case STATEVAR_PRIMITIVE: @@ -330,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/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 { diff --git a/examples/fluids/src/differential_filter.c b/examples/fluids/src/differential_filter.c index 7a05ff5d19..114dbc72b5 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; - - PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, num_qpts_1d, 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, &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); @@ -187,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)); @@ -198,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)); @@ -209,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 6cceabc156..2e19354629 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; @@ -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)); @@ -55,13 +56,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); - 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)); - 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, &basis_grid_aniso)); // -- Build RHS operator CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorProjection, AnisotropyTensorProjection_loc, &qf_rhs_assemble); @@ -136,17 +134,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); - - PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_qpts_1d, *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); @@ -156,7 +152,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/setupdm.c b/examples/fluids/src/setupdm.c index aa9912fb5d..8713f9eeb0 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" @@ -38,17 +39,45 @@ 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 + 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, q_order, &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, fe_coord_face_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 + 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)); + } + PetscCall(DMGetLabel(dm, "Face Sets", &label)); PetscCall(DMPlexLabelComplete(dm, label)); // Set wall BCs @@ -76,7 +105,13 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc 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)); + } + PetscCall(PetscFEDestroy(&fe)); } @@ -119,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, 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 0871d55234..4f763a1488 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -11,11 +11,9 @@ #include #include +#include #include "../navierstokes.h" -// 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) { @@ -38,14 +36,13 @@ 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) { DM dm_coord; - CeedInt Q_dim, loc_num_elem; + CeedInt loc_num_elem; PetscInt dim; CeedElemRestriction elem_restr_tmp; PetscFunctionBeginUser; PetscCall(DMGetDimension(dm, &dim)); dim -= height; - 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, 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, 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,6 +210,198 @@ 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(PETSC_SUCCESS); +} + +// ----------------------------------------------------------------------------- +// 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(PETSC_SUCCESS); +} + +// ----------------------------------------------------------------------------- +// Get CEED Basis from DMPlex +// ----------------------------------------------------------------------------- +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; + 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, CEED_GAUSS, basis); + } + + PetscFunctionReturn(PETSC_SUCCESS); +} + PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, AppCtx app_ctx, ProblemData *problem, SimpleBC bc) { PetscFunctionBeginUser; @@ -220,26 +409,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; // ----------------------------------------------------------------------------- // 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_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); // ----------------------------------------------------------------------------- // 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, 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, 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 +596,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_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)); // ----------------------------------------------------------------------------- // CEED QFunctions @@ -445,7 +642,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 26b136ece9..8cd11b281e 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; @@ -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,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, Q_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); @@ -151,8 +155,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; @@ -174,7 +177,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)); } } diff --git a/examples/fluids/src/turb_spanstats.c b/examples/fluids/src/turb_spanstats.c index 934d458d09..634226c615 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_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, &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, + 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); 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, &(*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, &(*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 7a84fd2763..75730096e1 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)); @@ -54,7 +55,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 +67,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); + PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, -1, 0, &elem_restr_grad_velo, NULL, NULL)); - PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, 0, num_qpts_1d, 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); + PetscCall(CreateBasisFromPlex(ceed, grad_velo_proj->dm, 0, 0, 0, 0, &basis_grad_velo)); // -- Build RHS operator switch (user->phys->state_var) { 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 de45762d35..5aee7101de 100644 Binary files a/examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin and b/examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin differ