1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 #include "../qfunctions/grid_anisotropy_tensor.h" 5 6 #include <petscdmplex.h> 7 8 #include <navierstokes.h> 9 10 PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso, 11 CeedVector *grid_aniso_vector) { 12 NodalProjectionData grid_aniso_proj; 13 OperatorApplyContext l2_rhs_ctx; 14 CeedOperator op_rhs_assemble, op_mass; 15 CeedQFunction qf_rhs_assemble, qf_mass; 16 CeedBasis basis_grid_aniso; 17 CeedInt q_data_size; 18 MPI_Comm comm = PetscObjectComm((PetscObject)user->dm); 19 KSP ksp; 20 DMLabel domain_label = NULL; 21 PetscInt label_value = 0, height = 0, dm_field = 0; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscNew(&grid_aniso_proj)); 25 26 // -- Create DM for Anisotropic tensor L^2 projection 27 grid_aniso_proj->num_comp = 7; 28 PetscCall(DMClone(user->dm, &grid_aniso_proj->dm)); 29 PetscCall(PetscObjectSetName((PetscObject)grid_aniso_proj->dm, "Grid Anisotropy Tensor Projection")); 30 31 { // -- Setup DM 32 PetscSection section; 33 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, 1, &grid_aniso_proj->num_comp, 34 grid_aniso_proj->dm)); 35 36 PetscCall(DMGetLocalSection(grid_aniso_proj->dm, §ion)); 37 PetscCall(PetscSectionSetFieldName(section, 0, "")); 38 PetscCall(PetscSectionSetComponentName(section, 0, 0, "KMGridAnisotropyTensorXX")); 39 PetscCall(PetscSectionSetComponentName(section, 0, 1, "KMGridAnisotropyTensorYY")); 40 PetscCall(PetscSectionSetComponentName(section, 0, 2, "KMGridAnisotropyTensorZZ")); 41 PetscCall(PetscSectionSetComponentName(section, 0, 3, "KMGridAnisotropyTensorYZ")); 42 PetscCall(PetscSectionSetComponentName(section, 0, 4, "KMGridAnisotropyTensorXZ")); 43 PetscCall(PetscSectionSetComponentName(section, 0, 5, "KMGridAnisotropyTensorXY")); 44 PetscCall(PetscSectionSetComponentName(section, 0, 6, "GridAnisotropyTensorFrobNorm")); 45 } 46 47 // -- Get Pre-requisite things 48 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size)); 49 50 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, grid_aniso_proj->dm, domain_label, label_value, height, dm_field, elem_restr_grid_aniso)); 51 PetscCall(CreateBasisFromPlex(ceed, grid_aniso_proj->dm, domain_label, label_value, height, dm_field, &basis_grid_aniso)); 52 53 // -- Build RHS operator 54 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorProjection, AnisotropyTensorProjection_loc, &qf_rhs_assemble)); 55 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "qdata", q_data_size, CEED_EVAL_NONE)); 56 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_assemble, "v", grid_aniso_proj->num_comp, CEED_EVAL_INTERP)); 57 58 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_assemble, NULL, NULL, &op_rhs_assemble)); 59 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 60 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE)); 61 62 PetscCall(OperatorApplyContextCreate(user->dm, grid_aniso_proj->dm, ceed, op_rhs_assemble, CEED_VECTOR_NONE, NULL, NULL, NULL, &l2_rhs_ctx)); 63 64 // -- Build Mass Operator 65 PetscCall(CreateMassQFunction(ceed, grid_aniso_proj->num_comp, q_data_size, &qf_mass)); 66 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 67 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE)); 68 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 69 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE)); 70 71 { // -- Setup KSP for L^2 projection 72 Mat mat_mass; 73 74 PetscCall(MatCreateCeed(grid_aniso_proj->dm, grid_aniso_proj->dm, op_mass, NULL, &mat_mass)); 75 76 PetscCall(KSPCreate(comm, &ksp)); 77 PetscCall(KSPSetOptionsPrefix(ksp, "grid_anisotropy_tensor_projection_")); 78 { 79 PC pc; 80 PetscCall(KSPGetPC(ksp, &pc)); 81 PetscCall(PCSetType(pc, PCJACOBI)); 82 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 83 PetscCall(KSPSetType(ksp, KSPCG)); 84 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 85 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 86 } 87 PetscCall(KSPSetFromOptions_WithMatCeed(ksp, mat_mass)); 88 PetscCall(MatDestroy(&mat_mass)); 89 } 90 91 { // -- Project anisotropy data and store in CeedVector 92 Vec Grid_Anisotropy, grid_anisotropy_loc; 93 94 // Get L^2 Projection RHS 95 PetscCall(DMGetGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 96 97 PetscCall(ApplyCeedOperatorLocalToGlobal(NULL, Grid_Anisotropy, l2_rhs_ctx)); 98 99 // Solve projection problem 100 PetscCall(KSPSolve(ksp, Grid_Anisotropy, Grid_Anisotropy)); 101 102 // Copy anisotropy tensor data to CeedVector 103 PetscCall(DMGetLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 104 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, grid_aniso_vector, NULL)); 105 PetscCall(DMGlobalToLocal(grid_aniso_proj->dm, Grid_Anisotropy, INSERT_VALUES, grid_anisotropy_loc)); 106 PetscCall(VecCopyPetscToCeed(grid_anisotropy_loc, *grid_aniso_vector)); 107 PetscCall(DMRestoreLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 108 PetscCall(DMRestoreGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 109 } 110 111 // -- Cleanup 112 PetscCall(NodalProjectionDataDestroy(grid_aniso_proj)); 113 PetscCall(OperatorApplyContextDestroy(l2_rhs_ctx)); 114 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_assemble)); 115 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 116 PetscCallCeed(ceed, CeedBasisDestroy(&basis_grid_aniso)); 117 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_assemble)); 118 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 119 PetscCall(KSPDestroy(&ksp)); 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso, 124 CeedVector *aniso_colloc_ceed, PetscInt *num_comp_aniso) { 125 CeedInt q_data_size, num_nodes; 126 CeedQFunction qf_colloc; 127 CeedOperator op_colloc; 128 DMLabel domain_label = NULL; 129 PetscInt label_value = 0, height = 0; 130 131 PetscFunctionBeginUser; 132 *num_comp_aniso = 7; 133 PetscCallCeed(ceed, CeedBasisGetNumNodes(ceed_data->basis_q, &num_nodes)); 134 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size)); 135 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, user->dm, domain_label, label_value, height, *num_comp_aniso, elem_restr_grid_aniso)); 136 137 // -- Build collocation operator 138 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc)); 139 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_colloc, "qdata", q_data_size, CEED_EVAL_NONE)); 140 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_colloc, "v", *num_comp_aniso, CEED_EVAL_NONE)); 141 142 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_colloc, NULL, NULL, &op_colloc)); 143 PetscCallCeed(ceed, CeedOperatorSetField(op_colloc, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 144 PetscCallCeed(ceed, CeedOperatorSetField(op_colloc, "v", *elem_restr_grid_aniso, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 145 146 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, aniso_colloc_ceed, NULL)); 147 148 PetscCallCeed(ceed, CeedOperatorApply(op_colloc, CEED_VECTOR_NONE, *aniso_colloc_ceed, CEED_REQUEST_IMMEDIATE)); 149 150 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_colloc)); 151 PetscCallCeed(ceed, CeedOperatorDestroy(&op_colloc)); 152 PetscFunctionReturn(PETSC_SUCCESS); 153 } 154