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/grid_anisotropy_tensor.h" 9 10 #include <petscdmplex.h> 11 12 #include "../navierstokes.h" 13 14 PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso, 15 CeedVector *grid_aniso_vector) { 16 NodalProjectionData grid_aniso_proj; 17 OperatorApplyContext mass_matop_ctx, l2_rhs_ctx; 18 CeedOperator op_rhs_assemble, op_mass; 19 CeedQFunction qf_rhs_assemble, qf_mass; 20 CeedBasis basis_grid_aniso; 21 PetscInt dim; 22 CeedInt num_qpts_1d, num_nodes_1d, q_data_size; 23 MPI_Comm comm = PetscObjectComm((PetscObject)user->dm); 24 KSP ksp; 25 26 PetscFunctionBeginUser; 27 PetscCall(PetscNew(&grid_aniso_proj)); 28 29 // -- Create DM for Anisotropic tensor L^2 projection 30 grid_aniso_proj->num_comp = 7; 31 PetscCall(DMClone(user->dm, &grid_aniso_proj->dm)); 32 PetscCall(DMGetDimension(grid_aniso_proj->dm, &dim)); 33 PetscCall(PetscObjectSetName((PetscObject)grid_aniso_proj->dm, "Grid Anisotropy Tensor Projection")); 34 35 { // -- Setup DM 36 PetscFE fe; 37 PetscSection section; 38 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grid_aniso_proj->num_comp, PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe)); 39 PetscCall(PetscObjectSetName((PetscObject)fe, "Grid Anisotropy Tensor Projection")); 40 PetscCall(DMAddField(grid_aniso_proj->dm, NULL, (PetscObject)fe)); 41 PetscCall(DMCreateDS(grid_aniso_proj->dm)); 42 PetscCall(DMPlexSetClosurePermutationTensor(grid_aniso_proj->dm, PETSC_DETERMINE, NULL)); 43 44 PetscCall(DMGetLocalSection(grid_aniso_proj->dm, §ion)); 45 PetscCall(PetscSectionSetFieldName(section, 0, "")); 46 PetscCall(PetscSectionSetComponentName(section, 0, 0, "KMGridAnisotropyTensorXX")); 47 PetscCall(PetscSectionSetComponentName(section, 0, 1, "KMGridAnisotropyTensorYY")); 48 PetscCall(PetscSectionSetComponentName(section, 0, 2, "KMGridAnisotropyTensorZZ")); 49 PetscCall(PetscSectionSetComponentName(section, 0, 3, "KMGridAnisotropyTensorYZ")); 50 PetscCall(PetscSectionSetComponentName(section, 0, 4, "KMGridAnisotropyTensorXZ")); 51 PetscCall(PetscSectionSetComponentName(section, 0, 5, "KMGridAnisotropyTensorXY")); 52 PetscCall(PetscSectionSetComponentName(section, 0, 6, "GridAnisotropyTensorFrobNorm")); 53 54 PetscCall(PetscFEDestroy(&fe)); 55 } 56 57 // -- Get Pre-requisite things 58 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); 59 CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); 60 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); 61 62 PetscCall( 63 GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, num_qpts_1d, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); 64 CeedBasisCreateTensorH1Lagrange(ceed, dim, grid_aniso_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grid_aniso); 65 66 // -- Build RHS operator 67 CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorProjection, AnisotropyTensorProjection_loc, &qf_rhs_assemble); 68 CeedQFunctionAddInput(qf_rhs_assemble, "qdata", q_data_size, CEED_EVAL_NONE); 69 CeedQFunctionAddOutput(qf_rhs_assemble, "v", grid_aniso_proj->num_comp, CEED_EVAL_INTERP); 70 71 CeedOperatorCreate(ceed, qf_rhs_assemble, NULL, NULL, &op_rhs_assemble); 72 CeedOperatorSetField(op_rhs_assemble, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 73 CeedOperatorSetField(op_rhs_assemble, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE); 74 75 PetscCall(OperatorApplyContextCreate(user->dm, grid_aniso_proj->dm, ceed, op_rhs_assemble, CEED_VECTOR_NONE, NULL, NULL, NULL, &l2_rhs_ctx)); 76 77 // -- Build Mass Operator 78 PetscCall(CreateMassQFunction(ceed, grid_aniso_proj->num_comp, q_data_size, &qf_mass)); 79 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 80 CeedOperatorSetField(op_mass, "u", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE); 81 CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 82 CeedOperatorSetField(op_mass, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE); 83 84 { // -- Setup KSP for L^2 projection 85 Mat mat_mass; 86 PetscCall(OperatorApplyContextCreate(grid_aniso_proj->dm, grid_aniso_proj->dm, ceed, op_mass, NULL, NULL, NULL, NULL, &mass_matop_ctx)); 87 PetscCall(CreateMatShell_Ceed(mass_matop_ctx, &mat_mass)); 88 89 PetscCall(KSPCreate(comm, &ksp)); 90 PetscCall(KSPSetOptionsPrefix(ksp, "grid_anisotropy_tensor_projection_")); 91 { 92 PC pc; 93 PetscCall(KSPGetPC(ksp, &pc)); 94 PetscCall(PCSetType(pc, PCJACOBI)); 95 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 96 PetscCall(KSPSetType(ksp, KSPCG)); 97 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 98 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 99 } 100 PetscCall(KSPSetOperators(ksp, mat_mass, mat_mass)); 101 PetscCall(KSPSetFromOptions(ksp)); 102 } 103 104 { // -- Project anisotropy data and store in CeedVector 105 Vec Grid_Anisotropy, grid_anisotropy_loc; 106 107 // Get L^2 Projection RHS 108 PetscCall(DMGetGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 109 110 PetscCall(ApplyCeedOperatorLocalToGlobal(NULL, Grid_Anisotropy, l2_rhs_ctx)); 111 112 // Solve projection problem 113 PetscCall(KSPSolve(ksp, Grid_Anisotropy, Grid_Anisotropy)); 114 115 // Copy anisotropy tensor data to CeedVector 116 PetscCall(DMGetLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 117 CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, grid_aniso_vector, NULL); 118 PetscCall(DMGlobalToLocal(grid_aniso_proj->dm, Grid_Anisotropy, INSERT_VALUES, grid_anisotropy_loc)); 119 PetscCall(VecCopyP2C(grid_anisotropy_loc, *grid_aniso_vector)); 120 PetscCall(DMRestoreLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 121 PetscCall(DMRestoreGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 122 } 123 124 // -- Cleanup 125 PetscCall(NodalProjectionDataDestroy(grid_aniso_proj)); 126 PetscCall(OperatorApplyContextDestroy(l2_rhs_ctx)); 127 CeedQFunctionDestroy(&qf_rhs_assemble); 128 CeedQFunctionDestroy(&qf_mass); 129 CeedBasisDestroy(&basis_grid_aniso); 130 CeedOperatorDestroy(&op_rhs_assemble); 131 CeedOperatorDestroy(&op_mass); 132 PetscCall(KSPDestroy(&ksp)); 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 136 PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso, 137 CeedVector *aniso_colloc_ceed, PetscInt *num_comp_aniso) { 138 PetscInt dim; 139 CeedInt q_data_size, num_qpts_1d, num_nodes_1d, loc_num_elem; 140 CeedQFunction qf_colloc; 141 CeedOperator op_colloc; 142 CeedBasis basis_grid_aniso; 143 144 PetscFunctionBeginUser; 145 // -- Get Pre-requisite things 146 *num_comp_aniso = 7; 147 PetscCall(DMGetDimension(user->dm, &dim)); 148 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); 149 CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); 150 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); 151 152 PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_qpts_1d, *num_comp_aniso, NULL, NULL, elem_restr_grid_aniso)); 153 154 CeedInt Q_dim = CeedIntPow(num_qpts_1d, dim); 155 CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &loc_num_elem); 156 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q_dim, *num_comp_aniso, *num_comp_aniso * loc_num_elem * Q_dim, CEED_STRIDES_BACKEND, 157 elem_restr_grid_aniso); 158 159 CeedBasisCreateTensorH1Lagrange(ceed, dim, *num_comp_aniso, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grid_aniso); 160 161 // -- Build collocation operator 162 CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc); 163 CeedQFunctionAddInput(qf_colloc, "qdata", q_data_size, CEED_EVAL_NONE); 164 CeedQFunctionAddOutput(qf_colloc, "v", *num_comp_aniso, CEED_EVAL_NONE); 165 166 CeedOperatorCreate(ceed, qf_colloc, NULL, NULL, &op_colloc); 167 CeedOperatorSetField(op_colloc, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 168 CeedOperatorSetField(op_colloc, "v", *elem_restr_grid_aniso, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 169 CeedOperatorSetNumQuadraturePoints(op_colloc, CeedIntPow(num_qpts_1d, dim)); 170 171 CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, aniso_colloc_ceed, NULL); 172 173 CeedOperatorApply(op_colloc, CEED_VECTOR_NONE, *aniso_colloc_ceed, CEED_REQUEST_IMMEDIATE); 174 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177