1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include "../qfunctions/sgs_dd_model.h" 9 10 #include <petscdmplex.h> 11 12 #include "../navierstokes.h" 13 14 typedef struct { 15 CeedElemRestriction elem_restr_grid_aniso, elem_restr_sgs; 16 CeedVector grid_aniso_ceed; 17 CeedQFunctionContext sgsdd_qfctx, ifunction_qfctx; 18 } *SgsDDSetupData; 19 20 PetscErrorCode SgsDDSetupDataDestroy(SgsDDSetupData sgs_dd_setup_data) { 21 Ceed ceed; 22 23 PetscFunctionBeginUser; 24 PetscCall(CeedElemRestrictionGetCeed(sgs_dd_setup_data->elem_restr_sgs, &ceed)); 25 26 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&sgs_dd_setup_data->elem_restr_grid_aniso)); 27 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&sgs_dd_setup_data->elem_restr_sgs)); 28 PetscCallCeed(ceed, CeedVectorDestroy(&sgs_dd_setup_data->grid_aniso_ceed)); 29 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&sgs_dd_setup_data->sgsdd_qfctx)); 30 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&sgs_dd_setup_data->ifunction_qfctx)); 31 PetscCall(PetscFree(sgs_dd_setup_data)); 32 PetscFunctionReturn(PETSC_SUCCESS); 33 } 34 35 // @brief Create DM for storing subgrid stress at nodes 36 static PetscErrorCode SgsDDCreateDM(DM dm_source, DM *dm_sgs, PetscInt degree, PetscInt q_extra, PetscInt *num_components) { 37 PetscSection section; 38 39 PetscFunctionBeginUser; 40 *num_components = 6; 41 42 PetscCall(DMClone(dm_source, dm_sgs)); 43 PetscCall(PetscObjectSetName((PetscObject)*dm_sgs, "Subgrid Stress Projection")); 44 45 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, num_components, *dm_sgs)); 46 47 PetscCall(DMGetLocalSection(*dm_sgs, §ion)); 48 PetscCall(PetscSectionSetFieldName(section, 0, "")); 49 PetscCall(PetscSectionSetComponentName(section, 0, 0, "KMSubgridStressXX")); 50 PetscCall(PetscSectionSetComponentName(section, 0, 1, "KMSubgridStressYY")); 51 PetscCall(PetscSectionSetComponentName(section, 0, 2, "KMSubgridStressZZ")); 52 PetscCall(PetscSectionSetComponentName(section, 0, 3, "KMSubgridStressYZ")); 53 PetscCall(PetscSectionSetComponentName(section, 0, 4, "KMSubgridStressXZ")); 54 PetscCall(PetscSectionSetComponentName(section, 0, 5, "KMSubgridStressXY")); 55 PetscFunctionReturn(PETSC_SUCCESS); 56 }; 57 58 // @brief Evaluate data-driven SGS using fused method 59 static PetscErrorCode SgsDDNodalStressEval_Fused(User user, Vec Q_loc, Vec VelocityGradient, Vec SGSNodal_loc) { 60 SgsDDData sgs_dd_data = user->sgs_dd_data; 61 PetscMemType q_mem_type; 62 63 PetscFunctionBeginUser; 64 PetscCall(VecPetscToCeed(Q_loc, &q_mem_type, user->q_ceed)); // q_ceed is an implicit input 65 66 PetscCall(ApplyCeedOperatorGlobalToLocal(VelocityGradient, SGSNodal_loc, sgs_dd_data->op_nodal_evaluation_ctx)); 67 68 PetscCall(VecCeedToPetsc(user->q_ceed, q_mem_type, Q_loc)); 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 72 // @brief Create CeedOperator to calculate data-drive SGS at nodes using fused operator 73 static PetscErrorCode SgsDDSetupNodalEvaluation_Fused(Ceed ceed, User user, CeedData ceed_data, SgsDDSetupData sgs_dd_setup_data) { 74 SgsDDData sgs_dd_data = user->sgs_dd_data; 75 CeedQFunction qf_sgs_dd_nodal; 76 CeedOperator op_sgs_dd_nodal; 77 CeedInt num_comp_q, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso; 78 PetscInt dim; 79 CeedVector inv_multiplicity; 80 CeedElemRestriction elem_restr_inv_multiplicity, elem_restr_grad_velo, elem_restr_sgs; 81 DMLabel domain_label = NULL; 82 PetscInt label_value = 0, height = 0, dm_field = 0; 83 84 PetscFunctionBeginUser; 85 PetscCall(DMGetDimension(user->dm, &dim)); 86 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x)); 87 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q)); 88 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(sgs_dd_setup_data->elem_restr_grid_aniso, &num_comp_grid_aniso)); 89 90 { // Get velocity gradient information 91 CeedOperatorField op_field; 92 PetscCallCeed(ceed, CeedOperatorGetFieldByName(user->grad_velo_proj->l2_rhs_ctx->op, "velocity gradient", &op_field)); 93 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_grad_velo)); 94 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_grad_velo, &num_comp_grad_velo)); 95 } 96 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, dm_field, &elem_restr_sgs)); 97 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_sgs, &sgs_dd_data->sgs_nodal_ceed, NULL)); 98 99 PetscCall(GetInverseMultiplicity(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, dm_field, PETSC_FALSE, &elem_restr_inv_multiplicity, 100 &inv_multiplicity)); 101 102 // -- Create operator for SGS DD model nodal evaluation 103 switch (user->phys->state_var) { 104 case STATEVAR_PRIMITIVE: 105 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDNodal_Prim, ComputeSgsDDNodal_Prim_loc, &qf_sgs_dd_nodal)); 106 break; 107 case STATEVAR_CONSERVATIVE: 108 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDNodal_Conserv, ComputeSgsDDNodal_Conserv_loc, &qf_sgs_dd_nodal)); 109 break; 110 default: 111 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Data-driven SGS nodal evaluation not available for chosen state variable"); 112 } 113 114 // Mesh/geometry order and solution basis order may differ, therefore must interpolate 115 CeedBasis basis_x_to_q; 116 PetscCallCeed(ceed, CeedBasisCreateProjection(ceed_data->basis_x, ceed_data->basis_q, &basis_x_to_q)); 117 118 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_sgs_dd_nodal, sgs_dd_setup_data->sgsdd_qfctx)); 119 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_nodal, "q", num_comp_q, CEED_EVAL_NONE)); 120 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_nodal, "x", num_comp_x, CEED_EVAL_INTERP)); 121 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_nodal, "gradient velocity", num_comp_grad_velo, CEED_EVAL_NONE)); 122 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_nodal, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE)); 123 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_nodal, "inverse multiplicity", 1, CEED_EVAL_NONE)); 124 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_dd_nodal, "km_sgs", sgs_dd_data->num_comp_sgs, CEED_EVAL_NONE)); 125 126 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_sgs_dd_nodal, NULL, NULL, &op_sgs_dd_nodal)); 127 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "q", ceed_data->elem_restr_q, CEED_BASIS_NONE, user->q_ceed)); 128 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "x", ceed_data->elem_restr_x, basis_x_to_q, ceed_data->x_coord)); 129 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "gradient velocity", elem_restr_grad_velo, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 130 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "anisotropy tensor", sgs_dd_setup_data->elem_restr_grid_aniso, CEED_BASIS_NONE, 131 sgs_dd_setup_data->grid_aniso_ceed)); 132 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "inverse multiplicity", elem_restr_inv_multiplicity, CEED_BASIS_NONE, inv_multiplicity)); 133 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "km_sgs", elem_restr_sgs, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 134 135 PetscCall(OperatorApplyContextCreate(user->grad_velo_proj->dm, sgs_dd_data->dm_sgs, ceed, op_sgs_dd_nodal, NULL, sgs_dd_data->sgs_nodal_ceed, NULL, 136 NULL, &sgs_dd_data->op_nodal_evaluation_ctx)); 137 138 sgs_dd_setup_data->elem_restr_sgs = elem_restr_sgs; 139 sgs_dd_data->sgs_nodal_eval = SgsDDNodalStressEval_Fused; 140 141 PetscCallCeed(ceed, CeedVectorDestroy(&inv_multiplicity)); 142 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q)); 143 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_inv_multiplicity)); 144 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_dd_nodal)); 145 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_dd_nodal)); 146 PetscFunctionReturn(PETSC_SUCCESS); 147 } 148 149 // @brief Setup data-driven model inference using internal (libCEED native) implementation 150 static PetscErrorCode SgsDDSetupNodalEvaluation_Sequential_Internal(Ceed ceed, SgsDDData sgs_dd_data, SgsDDSetupData sgs_dd_setup_data, 151 CeedElemRestriction elem_restr_dd_inputs, 152 CeedElemRestriction elem_restr_dd_outputs, 153 CeedElemRestriction elem_restr_inv_multiplicity, CeedVector inv_multiplicity, 154 void **ctx) { 155 CeedQFunction qf_sgs_dd_inference; 156 CeedOperator op_sgs_dd_inference; 157 OperatorApplyContext *op_context = (OperatorApplyContext *)ctx; 158 159 PetscFunctionBeginUser; 160 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDNodal_Sequential_Inference, ComputeSgsDDNodal_Sequential_Inference_loc, 161 &qf_sgs_dd_inference)); 162 163 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_sgs_dd_inference, sgs_dd_setup_data->sgsdd_qfctx)); 164 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_inference, "model inputs", sgs_dd_data->num_comp_inputs, CEED_EVAL_NONE)); 165 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_inference, "inverse multiplicity", 1, CEED_EVAL_NONE)); 166 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_dd_inference, "model outputs", sgs_dd_data->num_comp_outputs, CEED_EVAL_NONE)); 167 168 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_sgs_dd_inference, NULL, NULL, &op_sgs_dd_inference)); 169 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inference, "model inputs", elem_restr_dd_inputs, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 170 PetscCallCeed(ceed, 171 CeedOperatorSetField(op_sgs_dd_inference, "inverse multiplicity", elem_restr_inv_multiplicity, CEED_BASIS_NONE, inv_multiplicity)); 172 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inference, "model outputs", elem_restr_dd_outputs, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 173 174 PetscCall(OperatorApplyContextCreate(sgs_dd_data->dm_dd_inputs, sgs_dd_data->dm_dd_outputs, ceed, op_sgs_dd_inference, NULL, NULL, NULL, NULL, 175 op_context)); 176 sgs_dd_data->sgs_nodal_inference_ctx_destroy = (PetscErrorCode(*)(void *))OperatorApplyContextDestroy; 177 178 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_dd_inference)); 179 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_dd_inference)); 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 // @brief Perform data-driven model inference using internal (libCEED native) implementation 184 PetscErrorCode SgsDDNodalStressEval_Sequential_Internal(Vec DD_Inputs_loc, Vec DD_Outputs_loc, void *ctx) { 185 OperatorApplyContext op_context = *(OperatorApplyContext *)ctx; 186 187 PetscFunctionBeginUser; 188 PetscCall(ApplyCeedOperatorLocalToLocal(DD_Inputs_loc, DD_Outputs_loc, op_context)); 189 PetscFunctionReturn(PETSC_SUCCESS); 190 } 191 192 // @brief Evaluate data-driven SGS using sequential method 193 PetscErrorCode SgsDDNodalStressEval_Sequential(User user, Vec Q_loc, Vec VelocityGradient, Vec SGSNodal_loc) { 194 SgsDDData sgs_dd_data = user->sgs_dd_data; 195 PetscMemType q_mem_type; 196 Vec DD_Inputs_loc, DD_Outputs_loc; 197 198 PetscFunctionBeginUser; 199 PetscCall(DMGetLocalVector(sgs_dd_data->dm_dd_inputs, &DD_Inputs_loc)); 200 PetscCall(DMGetLocalVector(sgs_dd_data->dm_dd_outputs, &DD_Outputs_loc)); 201 PetscCall(VecPetscToCeed(Q_loc, &q_mem_type, user->q_ceed)); // q_ceed is an implicit input 202 203 PetscCall(ApplyCeedOperatorGlobalToLocal(VelocityGradient, DD_Inputs_loc, sgs_dd_data->op_nodal_dd_inputs_ctx)); 204 PetscCall(sgs_dd_data->sgs_nodal_inference(DD_Inputs_loc, DD_Outputs_loc, &sgs_dd_data->sgs_nodal_inference_ctx)); 205 PetscCall(ApplyCeedOperatorLocalToLocal(DD_Outputs_loc, SGSNodal_loc, sgs_dd_data->op_nodal_dd_outputs_ctx)); 206 207 PetscCall(VecCeedToPetsc(user->q_ceed, q_mem_type, Q_loc)); 208 PetscCall(DMRestoreLocalVector(sgs_dd_data->dm_dd_inputs, &DD_Inputs_loc)); 209 PetscCall(DMRestoreLocalVector(sgs_dd_data->dm_dd_outputs, &DD_Outputs_loc)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 // @brief Create CeedOperator to calculate data-drive SGS at nodes using sequentially-applied operators 214 static PetscErrorCode SgsDDSetupNodalEvaluation_Sequential(Ceed ceed, User user, CeedData ceed_data, SgsDDSetupData sgs_dd_setup_data) { 215 SgsDDData sgs_dd_data = user->sgs_dd_data; 216 CeedInt num_comp_q, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso, num_comp_eigvec = 9 + 1; 217 PetscInt dim; 218 CeedVector inv_multiplicity, eigvec; 219 CeedElemRestriction elem_restr_inv_multiplicity, elem_restr_grad_velo, elem_restr_sgs, elem_restr_eigvec, elem_restr_dd_inputs, 220 elem_restr_dd_outputs; 221 DMLabel domain_label = NULL; 222 PetscInt label_value = 0, height = 0, dm_field = 0; 223 224 PetscFunctionBeginUser; 225 { // Create DMs for data-driven input and output values 226 PetscSection section; 227 PetscInt degree, q_extra; 228 { // Get degree and number of quadrature points from dm_sgs 229 PetscFE fe; 230 PetscSpace basis; 231 PetscQuadrature quadrature; 232 PetscInt num_qpnts; 233 PetscCall(DMGetField(sgs_dd_data->dm_sgs, 0, NULL, (PetscObject *)&fe)); 234 PetscCall(PetscFEGetBasisSpace(fe, &basis)); 235 PetscCall(PetscSpaceGetDegree(basis, °ree, NULL)); 236 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 237 PetscCall(PetscQuadratureGetOrder(quadrature, &num_qpnts)); 238 q_extra = degree - num_qpnts; 239 } 240 241 PetscCall(DMClone(sgs_dd_data->dm_sgs, &sgs_dd_data->dm_dd_inputs)); 242 PetscCall(PetscObjectSetName((PetscObject)sgs_dd_data->dm_dd_inputs, "Data-Driven Model Inputs")); 243 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &sgs_dd_data->num_comp_inputs, sgs_dd_data->dm_dd_inputs)); 244 PetscCall(DMGetLocalSection(sgs_dd_data->dm_dd_inputs, §ion)); 245 PetscCall(PetscSectionSetFieldName(section, 0, "")); 246 for (CeedInt i = 0; i < sgs_dd_data->num_comp_inputs; i++) { 247 char component_name[PETSC_MAX_PATH_LEN]; 248 249 PetscCall(PetscSNPrintf(component_name, sizeof component_name, "DataDrivenInput%" CeedInt_FMT, i + 1)); 250 PetscCall(PetscSectionSetComponentName(section, 0, i, component_name)); 251 } 252 253 PetscCall(DMClone(sgs_dd_data->dm_sgs, &sgs_dd_data->dm_dd_outputs)); 254 PetscCall(PetscObjectSetName((PetscObject)sgs_dd_data->dm_dd_outputs, "Data-Driven Model Outputs")); 255 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &sgs_dd_data->num_comp_outputs, sgs_dd_data->dm_dd_outputs)); 256 PetscCall(DMGetLocalSection(sgs_dd_data->dm_dd_outputs, §ion)); 257 PetscCall(PetscSectionSetFieldName(section, 0, "")); 258 for (CeedInt i = 0; i < sgs_dd_data->num_comp_outputs; i++) { 259 char component_name[PETSC_MAX_PATH_LEN]; 260 261 PetscCall(PetscSNPrintf(component_name, sizeof component_name, "DataDrivenOutput%" CeedInt_FMT, i + 1)); 262 PetscCall(PetscSectionSetComponentName(section, 0, i, component_name)); 263 } 264 } 265 266 PetscCall(DMGetDimension(user->dm, &dim)); 267 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x)); 268 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q)); 269 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(sgs_dd_setup_data->elem_restr_grid_aniso, &num_comp_grid_aniso)); 270 271 { // Get velocity gradient information 272 CeedOperatorField op_field; 273 PetscCallCeed(ceed, CeedOperatorGetFieldByName(user->grad_velo_proj->l2_rhs_ctx->op, "velocity gradient", &op_field)); 274 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_grad_velo)); 275 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_grad_velo, &num_comp_grad_velo)); 276 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_grad_velo, &sgs_dd_data->grad_velo_ceed, NULL)); 277 } 278 279 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, dm_field, &elem_restr_sgs)); 280 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_sgs, &sgs_dd_data->sgs_nodal_ceed, NULL)); 281 PetscCall( 282 DMPlexCeedElemRestrictionCollocatedCreate(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, num_comp_eigvec, &elem_restr_eigvec)); 283 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_eigvec, &eigvec, NULL)); 284 285 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, sgs_dd_data->dm_dd_inputs, domain_label, label_value, height, dm_field, &elem_restr_dd_inputs)); 286 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, sgs_dd_data->dm_dd_outputs, domain_label, label_value, height, dm_field, &elem_restr_dd_outputs)); 287 288 PetscCall(GetInverseMultiplicity(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, dm_field, PETSC_FALSE, &elem_restr_inv_multiplicity, 289 &inv_multiplicity)); 290 291 { // Create operator for data-driven input evaluation 292 CeedQFunction qf_sgs_dd_inputs; 293 CeedOperator op_sgs_dd_inputs; 294 295 switch (user->phys->state_var) { 296 case STATEVAR_PRIMITIVE: 297 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDNodal_Sequential_Inputs_Prim, 298 ComputeSgsDDNodal_Sequential_Inputs_Prim_loc, &qf_sgs_dd_inputs)); 299 break; 300 case STATEVAR_CONSERVATIVE: 301 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDNodal_Sequential_Inputs_Conserv, 302 ComputeSgsDDNodal_Sequential_Inputs_Conserv_loc, &qf_sgs_dd_inputs)); 303 break; 304 default: 305 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, 306 "Data-driven SGS nodal input evaluation not available for chosen state variable"); 307 } 308 309 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_sgs_dd_inputs, sgs_dd_setup_data->sgsdd_qfctx)); 310 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_inputs, "q", num_comp_q, CEED_EVAL_NONE)); 311 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_inputs, "gradient velocity", num_comp_grad_velo, CEED_EVAL_NONE)); 312 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_inputs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE)); 313 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_inputs, "inverse multiplicity", 1, CEED_EVAL_NONE)); 314 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_dd_inputs, "eigenvectors", num_comp_eigvec, CEED_EVAL_NONE)); 315 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_dd_inputs, "model inputs", sgs_dd_data->num_comp_inputs, CEED_EVAL_NONE)); 316 317 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_sgs_dd_inputs, NULL, NULL, &op_sgs_dd_inputs)); 318 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inputs, "q", ceed_data->elem_restr_q, CEED_BASIS_NONE, user->q_ceed)); 319 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inputs, "gradient velocity", elem_restr_grad_velo, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 320 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inputs, "anisotropy tensor", sgs_dd_setup_data->elem_restr_grid_aniso, CEED_BASIS_NONE, 321 sgs_dd_setup_data->grid_aniso_ceed)); 322 PetscCallCeed(ceed, 323 CeedOperatorSetField(op_sgs_dd_inputs, "inverse multiplicity", elem_restr_inv_multiplicity, CEED_BASIS_NONE, inv_multiplicity)); 324 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inputs, "eigenvectors", elem_restr_eigvec, CEED_BASIS_NONE, eigvec)); 325 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_inputs, "model inputs", elem_restr_dd_inputs, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 326 327 PetscCall(OperatorApplyContextCreate(user->grad_velo_proj->dm, sgs_dd_data->dm_dd_inputs, ceed, op_sgs_dd_inputs, NULL, NULL, NULL, NULL, 328 &sgs_dd_data->op_nodal_dd_inputs_ctx)); 329 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_dd_inputs)); 330 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_dd_inputs)); 331 } 332 333 { // Create operator for data-driven output handling 334 CeedQFunction qf_sgs_dd_outputs; 335 CeedOperator op_sgs_dd_outputs; 336 337 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDNodal_Sequential_Outputs, ComputeSgsDDNodal_Sequential_Outputs_loc, 338 &qf_sgs_dd_outputs)); 339 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_sgs_dd_outputs, sgs_dd_setup_data->sgsdd_qfctx)); 340 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_outputs, "model outputs", sgs_dd_data->num_comp_outputs, CEED_EVAL_NONE)); 341 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_outputs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE)); 342 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_outputs, "inverse multiplicity", 1, CEED_EVAL_NONE)); 343 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_dd_outputs, "eigenvectors", num_comp_eigvec, CEED_EVAL_NONE)); 344 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_dd_outputs, "km_sgs", sgs_dd_data->num_comp_sgs, CEED_EVAL_NONE)); 345 346 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_sgs_dd_outputs, NULL, NULL, &op_sgs_dd_outputs)); 347 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_outputs, "model outputs", elem_restr_dd_outputs, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 348 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_outputs, "anisotropy tensor", sgs_dd_setup_data->elem_restr_grid_aniso, CEED_BASIS_NONE, 349 sgs_dd_setup_data->grid_aniso_ceed)); 350 PetscCallCeed(ceed, 351 CeedOperatorSetField(op_sgs_dd_outputs, "inverse multiplicity", elem_restr_inv_multiplicity, CEED_BASIS_NONE, inv_multiplicity)); 352 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_outputs, "eigenvectors", elem_restr_eigvec, CEED_BASIS_NONE, eigvec)); 353 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_outputs, "km_sgs", elem_restr_sgs, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 354 355 PetscCall(OperatorApplyContextCreate(sgs_dd_data->dm_dd_outputs, sgs_dd_data->dm_sgs, ceed, op_sgs_dd_outputs, NULL, sgs_dd_data->sgs_nodal_ceed, 356 NULL, NULL, &sgs_dd_data->op_nodal_dd_outputs_ctx)); 357 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_dd_outputs)); 358 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_dd_outputs)); 359 } 360 361 sgs_dd_data->sgs_nodal_inference = SgsDDNodalStressEval_Sequential_Internal; 362 sgs_dd_data->sgs_nodal_eval = SgsDDNodalStressEval_Sequential; 363 PetscCall(SgsDDSetupNodalEvaluation_Sequential_Internal(ceed, sgs_dd_data, sgs_dd_setup_data, elem_restr_dd_inputs, elem_restr_dd_outputs, 364 elem_restr_inv_multiplicity, inv_multiplicity, &sgs_dd_data->sgs_nodal_inference_ctx)); 365 366 sgs_dd_setup_data->elem_restr_sgs = elem_restr_sgs; 367 368 PetscCallCeed(ceed, CeedVectorDestroy(&inv_multiplicity)); 369 PetscCallCeed(ceed, CeedVectorDestroy(&eigvec)); 370 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_inv_multiplicity)); 371 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_eigvec)); 372 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_dd_inputs)); 373 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_dd_outputs)); 374 PetscFunctionReturn(PETSC_SUCCESS); 375 } 376 377 // @brief Create CeedOperator to compute SGS contribution to the residual 378 static PetscErrorCode SgsSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_data, SgsDDSetupData sgs_dd_setup_data) { 379 SgsDDData sgs_dd_data = user->sgs_dd_data; 380 CeedInt num_comp_q, num_comp_qd, num_comp_x; 381 PetscInt dim; 382 CeedQFunction qf_sgs_apply; 383 CeedOperator op_sgs_apply; 384 CeedBasis basis_sgs; 385 386 PetscFunctionBeginUser; 387 PetscCall(DMGetDimension(user->dm, &dim)); 388 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q)); 389 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd)); 390 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x)); 391 392 PetscCall(CreateBasisFromPlex(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, &basis_sgs)); 393 394 switch (user->phys->state_var) { 395 case STATEVAR_PRIMITIVE: 396 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, IFunction_NodalSgs_Prim, IFunction_NodalSgs_Prim_loc, &qf_sgs_apply)); 397 break; 398 case STATEVAR_CONSERVATIVE: 399 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, IFunction_NodalSgs_Conserv, IFunction_NodalSgs_Conserv_loc, &qf_sgs_apply)); 400 break; 401 default: 402 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Nodal SGS evaluation not available for chosen state variable"); 403 } 404 405 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_sgs_apply, sgs_dd_setup_data->ifunction_qfctx)); 406 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "q", num_comp_q, CEED_EVAL_INTERP)); 407 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "qdata", num_comp_qd, CEED_EVAL_NONE)); 408 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "km_sgs", sgs_dd_data->num_comp_sgs, CEED_EVAL_INTERP)); 409 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_apply, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD)); 410 411 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_sgs_apply, NULL, NULL, &op_sgs_apply)); 412 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 413 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 414 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "km_sgs", sgs_dd_setup_data->elem_restr_sgs, basis_sgs, sgs_dd_data->sgs_nodal_ceed)); 415 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "Grad_v", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 416 417 PetscCall( 418 OperatorApplyContextCreate(user->dm, user->dm, ceed, op_sgs_apply, user->q_ceed, user->g_ceed, NULL, NULL, &sgs_dd_data->op_sgs_apply_ctx)); 419 420 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_apply)); 421 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_apply)); 422 PetscFunctionReturn(PETSC_SUCCESS); 423 } 424 425 // @brief Calculate and add data-driven SGS residual to the global residual 426 PetscErrorCode SgsDDApplyIFunction(User user, const Vec Q_loc, Vec G_loc) { 427 SgsDDData sgs_dd_data = user->sgs_dd_data; 428 Vec VelocityGradient, SGSNodal_loc; 429 PetscMemType sgs_nodal_mem_type; 430 431 PetscFunctionBeginUser; 432 PetscCall(DMGetGlobalVector(user->grad_velo_proj->dm, &VelocityGradient)); 433 PetscCall(VelocityGradientProjectionApply(user->grad_velo_proj, Q_loc, VelocityGradient)); 434 435 // -- Compute Nodal SGS tensor 436 PetscCall(DMGetLocalVector(sgs_dd_data->dm_sgs, &SGSNodal_loc)); 437 PetscCall(sgs_dd_data->sgs_nodal_eval(user, Q_loc, VelocityGradient, SGSNodal_loc)); 438 439 // -- Compute contribution of the SGS stress 440 PetscCall(VecPetscToCeed(SGSNodal_loc, &sgs_nodal_mem_type, sgs_dd_data->sgs_nodal_ceed)); // sgs_nodal_ceed is an implicit input 441 PetscCall(ApplyAddCeedOperatorLocalToLocal(Q_loc, G_loc, sgs_dd_data->op_sgs_apply_ctx)); 442 443 // -- Return local SGS vector 444 PetscCall(VecCeedToPetsc(sgs_dd_data->sgs_nodal_ceed, sgs_nodal_mem_type, SGSNodal_loc)); 445 PetscCall(DMRestoreLocalVector(sgs_dd_data->dm_sgs, &SGSNodal_loc)); 446 PetscCall(DMRestoreGlobalVector(user->grad_velo_proj->dm, &VelocityGradient)); 447 PetscFunctionReturn(PETSC_SUCCESS); 448 } 449 450 // @brief B = A^T, A is NxM, B is MxN 451 static PetscErrorCode TransposeMatrix(const PetscScalar *A, PetscScalar *B, const PetscInt N, const PetscInt M) { 452 PetscFunctionBeginUser; 453 for (PetscInt i = 0; i < N; i++) { 454 for (PetscInt j = 0; j < M; j++) { 455 B[j * N + i] = A[i * M + j]; 456 } 457 } 458 PetscFunctionReturn(PETSC_SUCCESS); 459 } 460 461 // @brief Read neural network coefficients from file and put into context struct 462 static PetscErrorCode SgsDDContextFill(MPI_Comm comm, char data_dir[PETSC_MAX_PATH_LEN], SgsDDContext *psgsdd_ctx) { 463 SgsDDContext sgsdd_ctx; 464 PetscInt num_inputs = (*psgsdd_ctx)->num_inputs, num_outputs = (*psgsdd_ctx)->num_outputs, num_neurons = (*psgsdd_ctx)->num_neurons; 465 char file_path[PETSC_MAX_PATH_LEN]; 466 PetscScalar *temp; 467 468 PetscFunctionBeginUser; 469 { 470 SgsDDContext sgsdd_temp; 471 PetscCall(PetscNew(&sgsdd_temp)); 472 *sgsdd_temp = **psgsdd_ctx; 473 sgsdd_temp->offsets.bias1 = 0; 474 sgsdd_temp->offsets.bias2 = sgsdd_temp->offsets.bias1 + num_neurons; 475 sgsdd_temp->offsets.weight1 = sgsdd_temp->offsets.bias2 + num_neurons; 476 sgsdd_temp->offsets.weight2 = sgsdd_temp->offsets.weight1 + num_neurons * num_inputs; 477 sgsdd_temp->offsets.out_scaling = sgsdd_temp->offsets.weight2 + num_inputs * num_neurons; 478 PetscInt total_num_scalars = sgsdd_temp->offsets.out_scaling + 2 * num_outputs; 479 sgsdd_temp->total_bytes = sizeof(*sgsdd_ctx) + total_num_scalars * sizeof(sgsdd_ctx->data[0]); 480 PetscCall(PetscMalloc(sgsdd_temp->total_bytes, &sgsdd_ctx)); 481 *sgsdd_ctx = *sgsdd_temp; 482 PetscCall(PetscFree(sgsdd_temp)); 483 } 484 485 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "b1.dat")); 486 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, &sgsdd_ctx->data[sgsdd_ctx->offsets.bias1])); 487 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "b2.dat")); 488 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, &sgsdd_ctx->data[sgsdd_ctx->offsets.bias2])); 489 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "OutScaling.dat")); 490 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, &sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling])); 491 492 { 493 PetscCall(PetscMalloc1(num_inputs * num_neurons, &temp)); 494 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "w1.dat")); 495 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, temp)); 496 PetscCall(TransposeMatrix(temp, &sgsdd_ctx->data[sgsdd_ctx->offsets.weight1], num_inputs, num_neurons)); 497 PetscCall(PetscFree(temp)); 498 } 499 { 500 PetscCall(PetscMalloc1(num_outputs * num_neurons, &temp)); 501 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "w2.dat")); 502 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, temp)); 503 PetscCall(TransposeMatrix(temp, &sgsdd_ctx->data[sgsdd_ctx->offsets.weight2], num_neurons, num_outputs)); 504 PetscCall(PetscFree(temp)); 505 } 506 507 PetscCall(PetscFree(*psgsdd_ctx)); 508 *psgsdd_ctx = sgsdd_ctx; 509 PetscFunctionReturn(PETSC_SUCCESS); 510 } 511 512 PetscErrorCode SgsDDSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) { 513 PetscReal alpha = 0; 514 SgsDDContext sgsdd_ctx; 515 MPI_Comm comm = user->comm; 516 char sgs_dd_dir[PETSC_MAX_PATH_LEN] = "./dd_sgs_parameters"; 517 SgsDDSetupData sgs_dd_setup_data; 518 PetscBool use_fused; 519 NewtonianIdealGasContext gas; 520 521 PetscFunctionBeginUser; 522 PetscCall(VelocityGradientProjectionSetup(ceed, user, ceed_data, problem, user->phys->state_var, ceed_data->elem_restr_q, ceed_data->basis_q, 523 &user->grad_velo_proj)); 524 525 PetscCall(PetscNew(&user->sgs_dd_data)); 526 user->sgs_dd_data->num_comp_inputs = 6; 527 user->sgs_dd_data->num_comp_outputs = 6; 528 529 use_fused = PETSC_TRUE; 530 PetscOptionsBegin(comm, NULL, "SGS Data-Driven Model Options", NULL); 531 PetscCall(PetscOptionsReal("-sgs_model_dd_leakyrelu_alpha", "Slope parameter for Leaky ReLU activation function", NULL, alpha, &alpha, NULL)); 532 PetscCall(PetscOptionsString("-sgs_model_dd_parameter_dir", "Path to directory with model parameters (weights, biases, etc.)", NULL, sgs_dd_dir, 533 sgs_dd_dir, sizeof(sgs_dd_dir), NULL)); 534 PetscCall( 535 PetscOptionsBool("-sgs_model_dd_use_fused", "Use the fused SGS DD model evaluation instead of sequential", NULL, use_fused, &use_fused, NULL)); 536 PetscOptionsEnd(); 537 538 PetscCall(PetscNew(&sgsdd_ctx)); 539 sgsdd_ctx->num_layers = 1; 540 sgsdd_ctx->num_inputs = 6; 541 sgsdd_ctx->num_outputs = 6; 542 sgsdd_ctx->num_neurons = 20; 543 sgsdd_ctx->alpha = alpha; 544 545 PetscCall(SgsDDContextFill(comm, sgs_dd_dir, &sgsdd_ctx)); 546 547 // -- Create DM for storing SGS tensor at nodes 548 PetscCall(SgsDDCreateDM(user->dm, &user->sgs_dd_data->dm_sgs, user->app_ctx->degree, user->app_ctx->q_extra, &user->sgs_dd_data->num_comp_sgs)); 549 550 PetscCall(PetscNew(&sgs_dd_setup_data)); 551 552 PetscCallCeed(ceed, CeedQFunctionContextGetDataRead(problem->apply_vol_ifunction.qfunction_context, CEED_MEM_HOST, &gas)); 553 sgsdd_ctx->gas = *gas; 554 PetscCallCeed(ceed, CeedQFunctionContextRestoreDataRead(problem->apply_vol_ifunction.qfunction_context, &gas)); 555 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &sgs_dd_setup_data->sgsdd_qfctx)); 556 PetscCallCeed(ceed, 557 CeedQFunctionContextSetData(sgs_dd_setup_data->sgsdd_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sgsdd_ctx->total_bytes, sgsdd_ctx)); 558 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(sgs_dd_setup_data->sgsdd_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 559 560 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(problem->apply_vol_ifunction.qfunction_context, &sgs_dd_setup_data->ifunction_qfctx)); 561 562 // -- Compute and store anisotropy tensor 563 PetscCall(GridAnisotropyTensorProjectionSetupApply(ceed, user, ceed_data, &sgs_dd_setup_data->elem_restr_grid_aniso, 564 &sgs_dd_setup_data->grid_aniso_ceed)); 565 566 // -- Create Nodal Evaluation Operator 567 if (use_fused) PetscCall(SgsDDSetupNodalEvaluation_Fused(ceed, user, ceed_data, sgs_dd_setup_data)); 568 else PetscCall(SgsDDSetupNodalEvaluation_Sequential(ceed, user, ceed_data, sgs_dd_setup_data)); 569 570 // -- Create Operator to evalutate residual of SGS stress 571 PetscCall(SgsSetupNodalIFunction(ceed, user, ceed_data, sgs_dd_setup_data)); 572 573 PetscCall(SgsDDSetupDataDestroy(sgs_dd_setup_data)); 574 PetscFunctionReturn(PETSC_SUCCESS); 575 } 576 577 PetscErrorCode SgsDDDataDestroy(SgsDDData sgs_dd_data) { 578 PetscFunctionBeginUser; 579 if (!sgs_dd_data) PetscFunctionReturn(PETSC_SUCCESS); 580 Ceed ceed = sgs_dd_data->op_sgs_apply_ctx->ceed; 581 582 PetscCallCeed(ceed, CeedVectorDestroy(&sgs_dd_data->sgs_nodal_ceed)); 583 PetscCallCeed(ceed, CeedVectorDestroy(&sgs_dd_data->grad_velo_ceed)); 584 PetscCall(OperatorApplyContextDestroy(sgs_dd_data->op_nodal_evaluation_ctx)); 585 PetscCall(OperatorApplyContextDestroy(sgs_dd_data->op_sgs_apply_ctx)); 586 PetscCall(OperatorApplyContextDestroy(sgs_dd_data->op_nodal_dd_inputs_ctx)); 587 PetscCall(OperatorApplyContextDestroy(sgs_dd_data->op_nodal_dd_outputs_ctx)); 588 PetscCall(DMDestroy(&sgs_dd_data->dm_sgs)); 589 PetscCall(DMDestroy(&sgs_dd_data->dm_dd_inputs)); 590 PetscCall(DMDestroy(&sgs_dd_data->dm_dd_outputs)); 591 if (sgs_dd_data->sgs_nodal_inference_ctx) PetscCall(sgs_dd_data->sgs_nodal_inference_ctx_destroy(sgs_dd_data->sgs_nodal_inference_ctx)); 592 PetscCall(PetscFree(sgs_dd_data)); 593 PetscFunctionReturn(PETSC_SUCCESS); 594 } 595