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