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