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