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