1 // Copyright (c) 2017-2023, 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; 18 } *SgsDDModelSetupData; 19 20 PetscErrorCode SgsDDModelSetupDataDestroy(SgsDDModelSetupData sgs_dd_setup_data) { 21 Ceed ceed; 22 PetscFunctionBeginUser; 23 PetscCall(CeedElemRestrictionGetCeed(sgs_dd_setup_data->elem_restr_sgs, &ceed)); 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 29 PetscCall(PetscFree(sgs_dd_setup_data)); 30 PetscFunctionReturn(PETSC_SUCCESS); 31 } 32 33 // @brief Create DM for storing subgrid stress at nodes 34 PetscErrorCode SgsDDModelCreateDM(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(PetscObjectSetName((PetscObject)*dm_sgs, "Subgrid Stress Projection")); 42 43 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, num_components, *dm_sgs)); 44 45 PetscCall(DMGetLocalSection(*dm_sgs, §ion)); 46 PetscCall(PetscSectionSetFieldName(section, 0, "")); 47 PetscCall(PetscSectionSetComponentName(section, 0, 0, "KMSubgridStressXX")); 48 PetscCall(PetscSectionSetComponentName(section, 0, 1, "KMSubgridStressYY")); 49 PetscCall(PetscSectionSetComponentName(section, 0, 2, "KMSubgridStressZZ")); 50 PetscCall(PetscSectionSetComponentName(section, 0, 3, "KMSubgridStressYZ")); 51 PetscCall(PetscSectionSetComponentName(section, 0, 4, "KMSubgridStressXZ")); 52 PetscCall(PetscSectionSetComponentName(section, 0, 5, "KMSubgridStressXY")); 53 PetscFunctionReturn(PETSC_SUCCESS); 54 }; 55 56 // @brief Create CeedOperator to calculate data-drive SGS at nodes 57 PetscErrorCode SgsDDModelSetupNodalEvaluation(Ceed ceed, User user, CeedData ceed_data, SgsDDModelSetupData sgs_dd_setup_data) { 58 SgsDDData sgs_dd_data = user->sgs_dd_data; 59 CeedQFunction qf_multiplicity, qf_sgs_dd_nodal; 60 CeedOperator op_multiplicity, op_sgs_dd_nodal; 61 CeedInt num_comp_q, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso; 62 PetscInt dim; 63 CeedVector multiplicity, inv_multiplicity; 64 CeedElemRestriction elem_restr_inv_multiplicity, elem_restr_grad_velo, elem_restr_sgs; 65 DMLabel domain_label = NULL; 66 PetscInt label_value = 0, height = 0, dm_field = 0; 67 68 PetscFunctionBeginUser; 69 PetscCall(DMGetDimension(user->dm, &dim)); 70 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x)); 71 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q)); 72 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(sgs_dd_setup_data->elem_restr_grid_aniso, &num_comp_grid_aniso)); 73 74 { // Get velocity gradient information 75 CeedOperatorField op_field; 76 PetscCallCeed(ceed, CeedOperatorGetFieldByName(user->grad_velo_proj->l2_rhs_ctx->op, "velocity gradient", &op_field)); 77 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_grad_velo)); 78 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_grad_velo, &num_comp_grad_velo)); 79 } 80 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, dm_field, &elem_restr_sgs)); 81 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_sgs, &sgs_dd_data->sgs_nodal_ceed, NULL)); 82 83 // -- Create inverse multiplicity for correcting nodal assembly 84 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &multiplicity, NULL)); 85 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(ceed_data->elem_restr_q, multiplicity)); 86 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, sgs_dd_data->dm_sgs, domain_label, label_value, height, 1, &elem_restr_inv_multiplicity)); 87 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_inv_multiplicity, &inv_multiplicity, NULL)); 88 89 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, InverseMultiplicity, InverseMultiplicity_loc, &qf_multiplicity)); 90 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_multiplicity, "multiplicity", num_comp_q, CEED_EVAL_NONE)); 91 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_multiplicity, "inverse multiplicity", 1, CEED_EVAL_NONE)); 92 93 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_multiplicity, NULL, NULL, &op_multiplicity)); 94 PetscCallCeed(ceed, CeedOperatorSetName(op_multiplicity, "SGS DD Model - Create Multiplicity Scaling")); 95 PetscCallCeed(ceed, CeedOperatorSetField(op_multiplicity, "multiplicity", ceed_data->elem_restr_q, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE)); 96 PetscCallCeed( 97 ceed, CeedOperatorSetField(op_multiplicity, "inverse multiplicity", elem_restr_inv_multiplicity, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE)); 98 99 PetscCallCeed(ceed, CeedOperatorApply(op_multiplicity, multiplicity, inv_multiplicity, CEED_REQUEST_IMMEDIATE)); 100 101 // -- Create operator for SGS DD model nodal evaluation 102 switch (user->phys->state_var) { 103 case STATEVAR_PRIMITIVE: 104 PetscCallCeed(ceed, 105 CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDAnisotropicNodal_Prim, ComputeSgsDDAnisotropicNodal_Prim_loc, &qf_sgs_dd_nodal)); 106 break; 107 case STATEVAR_CONSERVATIVE: 108 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeSgsDDAnisotropicNodal_Conserv, ComputeSgsDDAnisotropicNodal_Conserv_loc, 109 &qf_sgs_dd_nodal)); 110 break; 111 default: 112 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, 113 "Anisotropic 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_COLLOCATED, 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_COLLOCATED, CEED_VECTOR_ACTIVE)); 132 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "anisotropy tensor", sgs_dd_setup_data->elem_restr_grid_aniso, CEED_BASIS_COLLOCATED, 133 sgs_dd_setup_data->grid_aniso_ceed)); 134 PetscCallCeed(ceed, 135 CeedOperatorSetField(op_sgs_dd_nodal, "inverse multiplicity", elem_restr_inv_multiplicity, CEED_BASIS_COLLOCATED, inv_multiplicity)); 136 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_dd_nodal, "km_sgs", elem_restr_sgs, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE)); 137 138 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, 139 NULL, &sgs_dd_data->op_nodal_evaluation_ctx)); 140 141 sgs_dd_setup_data->elem_restr_sgs = elem_restr_sgs; 142 143 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity)); 144 PetscCallCeed(ceed, CeedVectorDestroy(&inv_multiplicity)); 145 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q)); 146 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_inv_multiplicity)); 147 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_multiplicity)); 148 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_dd_nodal)); 149 PetscCallCeed(ceed, CeedOperatorDestroy(&op_multiplicity)); 150 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_dd_nodal)); 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 // @brief Create CeedOperator to compute SGS contribution to the residual 155 PetscErrorCode SgsModelSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_data, SgsDDModelSetupData sgs_dd_setup_data) { 156 SgsDDData sgs_dd_data = user->sgs_dd_data; 157 CeedInt num_comp_q, num_comp_qd, num_comp_x; 158 PetscInt dim; 159 CeedQFunction qf_sgs_apply; 160 CeedOperator op_sgs_apply; 161 CeedBasis basis_sgs; 162 163 PetscFunctionBeginUser; 164 PetscCall(DMGetDimension(user->dm, &dim)); 165 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q)); 166 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd)); 167 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x)); 168 169 PetscCall(CreateBasisFromPlex(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, &basis_sgs)); 170 171 switch (user->phys->state_var) { 172 case STATEVAR_PRIMITIVE: 173 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, IFunction_NodalSgs_Prim, IFunction_NodalSgs_Prim_loc, &qf_sgs_apply)); 174 break; 175 case STATEVAR_CONSERVATIVE: 176 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, IFunction_NodalSgs_Conserv, IFunction_NodalSgs_Conserv_loc, &qf_sgs_apply)); 177 break; 178 default: 179 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Nodal SGS evaluation not available for chosen state variable"); 180 } 181 182 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_sgs_apply, sgs_dd_setup_data->sgsdd_qfctx)); 183 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "q", num_comp_q, CEED_EVAL_INTERP)); 184 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "qdata", num_comp_qd, CEED_EVAL_NONE)); 185 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "x", num_comp_x, CEED_EVAL_INTERP)); 186 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_sgs_apply, "km_sgs", sgs_dd_data->num_comp_sgs, CEED_EVAL_INTERP)); 187 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_sgs_apply, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD)); 188 189 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_sgs_apply, NULL, NULL, &op_sgs_apply)); 190 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 191 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data)); 192 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord)); 193 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "km_sgs", sgs_dd_setup_data->elem_restr_sgs, basis_sgs, sgs_dd_data->sgs_nodal_ceed)); 194 PetscCallCeed(ceed, CeedOperatorSetField(op_sgs_apply, "Grad_v", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 195 196 PetscCall( 197 OperatorApplyContextCreate(user->dm, user->dm, ceed, op_sgs_apply, user->q_ceed, user->g_ceed, NULL, NULL, &sgs_dd_data->op_sgs_apply_ctx)); 198 199 PetscCallCeed(ceed, CeedOperatorDestroy(&op_sgs_apply)); 200 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_sgs_apply)); 201 PetscFunctionReturn(PETSC_SUCCESS); 202 } 203 204 // @brief Calculate and add data-driven SGS residual to the global residual 205 PetscErrorCode SgsDDModelApplyIFunction(User user, const Vec Q_loc, Vec G_loc) { 206 SgsDDData sgs_dd_data = user->sgs_dd_data; 207 Vec VelocityGradient, SGSNodal_loc; 208 PetscMemType sgs_nodal_mem_type, q_mem_type; 209 210 PetscFunctionBeginUser; 211 PetscCall(DMGetGlobalVector(user->grad_velo_proj->dm, &VelocityGradient)); 212 PetscCall(VelocityGradientProjectionApply(user, Q_loc, VelocityGradient)); 213 214 // -- Compute Nodal SGS tensor 215 PetscCall(DMGetLocalVector(sgs_dd_data->dm_sgs, &SGSNodal_loc)); 216 PetscCall(VecP2C(Q_loc, &q_mem_type, user->q_ceed)); // q_ceed is an implicit input 217 218 PetscCall(ApplyCeedOperatorGlobalToLocal(VelocityGradient, SGSNodal_loc, sgs_dd_data->op_nodal_evaluation_ctx)); 219 220 PetscCall(VecC2P(user->q_ceed, q_mem_type, Q_loc)); 221 PetscCall(VecP2C(SGSNodal_loc, &sgs_nodal_mem_type, sgs_dd_data->sgs_nodal_ceed)); // sgs_nodal_ceed is an implicit input 222 223 // -- Compute contribution of the SGS stress 224 PetscCall(ApplyAddCeedOperatorLocalToLocal(Q_loc, G_loc, sgs_dd_data->op_sgs_apply_ctx)); 225 226 // -- Return local SGS vector 227 PetscCall(VecC2P(sgs_dd_data->sgs_nodal_ceed, sgs_nodal_mem_type, SGSNodal_loc)); 228 PetscCall(DMRestoreLocalVector(sgs_dd_data->dm_sgs, &SGSNodal_loc)); 229 PetscCall(DMRestoreGlobalVector(user->grad_velo_proj->dm, &VelocityGradient)); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 // @brief B = A^T, A is NxM, B is MxN 234 PetscErrorCode TransposeMatrix(const PetscScalar *A, PetscScalar *B, const PetscInt N, const PetscInt M) { 235 PetscFunctionBeginUser; 236 for (PetscInt i = 0; i < N; i++) { 237 for (PetscInt j = 0; j < M; j++) { 238 B[j * N + i] = A[i * M + j]; 239 } 240 } 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 // @brief Read neural network coefficients from file and put into context struct 245 PetscErrorCode SgsDDModelContextFill(MPI_Comm comm, char data_dir[PETSC_MAX_PATH_LEN], SgsDDModelContext *psgsdd_ctx) { 246 SgsDDModelContext sgsdd_ctx; 247 PetscInt num_inputs = (*psgsdd_ctx)->num_inputs, num_outputs = (*psgsdd_ctx)->num_outputs, num_neurons = (*psgsdd_ctx)->num_neurons; 248 char file_path[PETSC_MAX_PATH_LEN]; 249 PetscScalar *temp; 250 251 PetscFunctionBeginUser; 252 { 253 SgsDDModelContext sgsdd_temp; 254 PetscCall(PetscNew(&sgsdd_temp)); 255 *sgsdd_temp = **psgsdd_ctx; 256 sgsdd_temp->offsets.bias1 = 0; 257 sgsdd_temp->offsets.bias2 = sgsdd_temp->offsets.bias1 + num_neurons; 258 sgsdd_temp->offsets.weight1 = sgsdd_temp->offsets.bias2 + num_neurons; 259 sgsdd_temp->offsets.weight2 = sgsdd_temp->offsets.weight1 + num_neurons * num_inputs; 260 sgsdd_temp->offsets.out_scaling = sgsdd_temp->offsets.weight2 + num_inputs * num_neurons; 261 PetscInt total_num_scalars = sgsdd_temp->offsets.out_scaling + 2 * num_outputs; 262 sgsdd_temp->total_bytes = sizeof(*sgsdd_ctx) + total_num_scalars * sizeof(sgsdd_ctx->data[0]); 263 PetscCall(PetscMalloc(sgsdd_temp->total_bytes, &sgsdd_ctx)); 264 *sgsdd_ctx = *sgsdd_temp; 265 PetscCall(PetscFree(sgsdd_temp)); 266 } 267 268 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "b1.dat")); 269 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, &sgsdd_ctx->data[sgsdd_ctx->offsets.bias1])); 270 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "b2.dat")); 271 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, &sgsdd_ctx->data[sgsdd_ctx->offsets.bias2])); 272 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "OutScaling.dat")); 273 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, &sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling])); 274 275 { 276 PetscCall(PetscMalloc1(num_inputs * num_neurons, &temp)); 277 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "w1.dat")); 278 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, temp)); 279 PetscCall(TransposeMatrix(temp, &sgsdd_ctx->data[sgsdd_ctx->offsets.weight1], num_inputs, num_neurons)); 280 PetscCall(PetscFree(temp)); 281 } 282 { 283 PetscCall(PetscMalloc1(num_outputs * num_neurons, &temp)); 284 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/%s", data_dir, "w2.dat")); 285 PetscCall(PhastaDatFileReadToArrayReal(comm, file_path, temp)); 286 PetscCall(TransposeMatrix(temp, &sgsdd_ctx->data[sgsdd_ctx->offsets.weight2], num_neurons, num_outputs)); 287 PetscCall(PetscFree(temp)); 288 } 289 290 PetscCall(PetscFree(*psgsdd_ctx)); 291 *psgsdd_ctx = sgsdd_ctx; 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 295 PetscErrorCode SgsDDModelSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 296 PetscReal alpha = 0; 297 SgsDDModelContext sgsdd_ctx; 298 MPI_Comm comm = user->comm; 299 char sgs_dd_dir[PETSC_MAX_PATH_LEN] = "./dd_sgs_parameters"; 300 SgsDDModelSetupData sgs_dd_setup_data; 301 NewtonianIdealGasContext gas; 302 303 PetscFunctionBeginUser; 304 PetscCall(VelocityGradientProjectionSetup(ceed, user, ceed_data, problem)); 305 306 PetscCall(PetscNew(&sgsdd_ctx)); 307 308 PetscOptionsBegin(comm, NULL, "SGS Data-Driven Model Options", NULL); 309 PetscCall(PetscOptionsReal("-sgs_model_dd_leakyrelu_alpha", "Slope parameter for Leaky ReLU activation function", NULL, alpha, &alpha, NULL)); 310 PetscCall(PetscOptionsString("-sgs_model_dd_parameter_dir", "Path to directory with model parameters (weights, biases, etc.)", NULL, sgs_dd_dir, 311 sgs_dd_dir, sizeof(sgs_dd_dir), NULL)); 312 PetscOptionsEnd(); 313 314 sgsdd_ctx->num_layers = 1; 315 sgsdd_ctx->num_inputs = 6; 316 sgsdd_ctx->num_outputs = 6; 317 sgsdd_ctx->num_neurons = 20; 318 sgsdd_ctx->alpha = alpha; 319 320 PetscCall(SgsDDModelContextFill(comm, sgs_dd_dir, &sgsdd_ctx)); 321 322 // -- Create DM for storing SGS tensor at nodes 323 PetscCall(PetscNew(&user->sgs_dd_data)); 324 PetscCall( 325 SgsDDModelCreateDM(user->dm, &user->sgs_dd_data->dm_sgs, user->app_ctx->degree, user->app_ctx->q_extra, &user->sgs_dd_data->num_comp_sgs)); 326 327 PetscCall(PetscNew(&sgs_dd_setup_data)); 328 329 PetscCallCeed(ceed, CeedQFunctionContextGetDataRead(problem->apply_vol_ifunction.qfunction_context, CEED_MEM_HOST, &gas)); 330 sgsdd_ctx->gas = *gas; 331 PetscCallCeed(ceed, CeedQFunctionContextRestoreDataRead(problem->apply_vol_ifunction.qfunction_context, &gas)); 332 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &sgs_dd_setup_data->sgsdd_qfctx)); 333 PetscCallCeed(ceed, 334 CeedQFunctionContextSetData(sgs_dd_setup_data->sgsdd_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sgsdd_ctx->total_bytes, sgsdd_ctx)); 335 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(sgs_dd_setup_data->sgsdd_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 336 337 // -- Compute and store anisotropy tensor 338 PetscCall(GridAnisotropyTensorProjectionSetupApply(ceed, user, ceed_data, &sgs_dd_setup_data->elem_restr_grid_aniso, 339 &sgs_dd_setup_data->grid_aniso_ceed)); 340 341 // -- Create Nodal Evaluation Operator 342 PetscCall(SgsDDModelSetupNodalEvaluation(ceed, user, ceed_data, sgs_dd_setup_data)); 343 344 // -- Create Operator to evalutate residual of SGS stress 345 PetscCall(SgsModelSetupNodalIFunction(ceed, user, ceed_data, sgs_dd_setup_data)); 346 347 PetscCall(SgsDDModelSetupDataDestroy(sgs_dd_setup_data)); 348 PetscFunctionReturn(PETSC_SUCCESS); 349 } 350 351 PetscErrorCode SgsDDDataDestroy(SgsDDData sgs_dd_data) { 352 PetscFunctionBeginUser; 353 if (!sgs_dd_data) PetscFunctionReturn(PETSC_SUCCESS); 354 Ceed ceed = sgs_dd_data->op_sgs_apply_ctx->ceed; 355 356 PetscCallCeed(ceed, CeedVectorDestroy(&sgs_dd_data->sgs_nodal_ceed)); 357 PetscCall(OperatorApplyContextDestroy(sgs_dd_data->op_nodal_evaluation_ctx)); 358 PetscCall(DMDestroy(&sgs_dd_data->dm_sgs)); 359 PetscCall(PetscFree(sgs_dd_data)); 360 PetscFunctionReturn(PETSC_SUCCESS); 361 } 362