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