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 CeedVector rhs_ceed; 22 PetscInt dim, q_data_size, num_qpts_1d, num_nodes_1d; 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(GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 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 CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, &rhs_ceed, NULL); 75 76 PetscCall(OperatorApplyContextCreate(user->dm, grid_aniso_proj->dm, ceed, op_rhs_assemble, ceed_data->q_data, rhs_ceed, NULL, NULL, &l2_rhs_ctx)); 77 78 // -- Build Mass Operator 79 PetscCall(CreateMassQFunction(ceed, grid_aniso_proj->num_comp, q_data_size, &qf_mass)); 80 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 81 CeedOperatorSetField(op_mass, "u", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE); 82 CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 83 CeedOperatorSetField(op_mass, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE); 84 85 { // -- Setup KSP for L^2 projection 86 PetscInt l_size, g_size; 87 Mat mat_mass; 88 VecType vec_type; 89 Vec M_inv; 90 CeedVector mass_output; 91 92 PetscCall(DMGetGlobalVector(grid_aniso_proj->dm, &M_inv)); 93 PetscCall(VecGetLocalSize(M_inv, &l_size)); 94 PetscCall(VecGetSize(M_inv, &g_size)); 95 PetscCall(VecGetType(M_inv, &vec_type)); 96 PetscCall(DMRestoreGlobalVector(grid_aniso_proj->dm, &M_inv)); 97 98 CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, &mass_output, NULL); 99 PetscCall( 100 OperatorApplyContextCreate(grid_aniso_proj->dm, grid_aniso_proj->dm, ceed, op_mass, rhs_ceed, mass_output, NULL, NULL, &mass_matop_ctx)); 101 CeedVectorDestroy(&mass_output); 102 103 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, mass_matop_ctx, &mat_mass)); 104 PetscCall(MatShellSetContextDestroy(mat_mass, (PetscErrorCode(*)(void *))OperatorApplyContextDestroy)); 105 PetscCall(MatShellSetOperation(mat_mass, MATOP_MULT, (void (*)(void))MatMult_Ceed)); 106 PetscCall(MatShellSetOperation(mat_mass, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiag_Ceed)); 107 PetscCall(MatShellSetVecType(mat_mass, vec_type)); 108 109 PetscCall(KSPCreate(comm, &ksp)); 110 PetscCall(KSPSetOptionsPrefix(ksp, "grid_anisotropy_tensor_projection_")); 111 { 112 PC pc; 113 PetscCall(KSPGetPC(ksp, &pc)); 114 PetscCall(PCSetType(pc, PCJACOBI)); 115 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 116 PetscCall(KSPSetType(ksp, KSPCG)); 117 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 118 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 119 } 120 PetscCall(KSPSetOperators(ksp, mat_mass, mat_mass)); 121 PetscCall(KSPSetFromOptions(ksp)); 122 } 123 124 { // -- Project anisotropy data and store in CeedVector 125 Vec Grid_Anisotropy, grid_anisotropy_loc; 126 PetscMemType mem_type; 127 128 // Get L^2 Projection RHS 129 PetscCall(DMGetGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 130 PetscCall(DMGetLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 131 132 PetscCall(VecP2C(grid_anisotropy_loc, &mem_type, l2_rhs_ctx->y_ceed)); 133 CeedOperatorApply(l2_rhs_ctx->op, CEED_VECTOR_NONE, l2_rhs_ctx->y_ceed, CEED_REQUEST_IMMEDIATE); 134 PetscCall(VecC2P(l2_rhs_ctx->y_ceed, mem_type, grid_anisotropy_loc)); 135 PetscCall(DMLocalToGlobal(l2_rhs_ctx->dm_y, grid_anisotropy_loc, ADD_VALUES, Grid_Anisotropy)); 136 137 // Solve projection problem 138 PetscCall(KSPSolve(ksp, Grid_Anisotropy, Grid_Anisotropy)); 139 140 // Copy anisotropy tensor data to CeedVector 141 CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, grid_aniso_vector, NULL); 142 PetscCall(DMGlobalToLocal(grid_aniso_proj->dm, Grid_Anisotropy, INSERT_VALUES, grid_anisotropy_loc)); 143 PetscCall(VecCopyP2C(grid_anisotropy_loc, *grid_aniso_vector)); 144 PetscCall(DMRestoreLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 145 PetscCall(DMRestoreGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 146 } 147 148 // -- Cleanup 149 PetscCall(NodalProjectionDataDestroy(grid_aniso_proj)); 150 PetscCall(OperatorApplyContextDestroy(l2_rhs_ctx)); 151 CeedVectorDestroy(&rhs_ceed); 152 CeedQFunctionDestroy(&qf_rhs_assemble); 153 CeedQFunctionDestroy(&qf_mass); 154 CeedBasisDestroy(&basis_grid_aniso); 155 CeedOperatorDestroy(&op_rhs_assemble); 156 CeedOperatorDestroy(&op_mass); 157 PetscCall(KSPDestroy(&ksp)); 158 PetscFunctionReturn(0); 159 } 160