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