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 CeedInt q_data_size; 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(PetscObjectSetName((PetscObject)grid_aniso_proj->dm, "Grid Anisotropy Tensor Projection")); 32 33 { // -- Setup DM 34 PetscSection section; 35 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, 1, &grid_aniso_proj->num_comp, 36 grid_aniso_proj->dm)); 37 38 PetscCall(DMGetLocalSection(grid_aniso_proj->dm, §ion)); 39 PetscCall(PetscSectionSetFieldName(section, 0, "")); 40 PetscCall(PetscSectionSetComponentName(section, 0, 0, "KMGridAnisotropyTensorXX")); 41 PetscCall(PetscSectionSetComponentName(section, 0, 1, "KMGridAnisotropyTensorYY")); 42 PetscCall(PetscSectionSetComponentName(section, 0, 2, "KMGridAnisotropyTensorZZ")); 43 PetscCall(PetscSectionSetComponentName(section, 0, 3, "KMGridAnisotropyTensorYZ")); 44 PetscCall(PetscSectionSetComponentName(section, 0, 4, "KMGridAnisotropyTensorXZ")); 45 PetscCall(PetscSectionSetComponentName(section, 0, 5, "KMGridAnisotropyTensorXY")); 46 PetscCall(PetscSectionSetComponentName(section, 0, 6, "GridAnisotropyTensorFrobNorm")); 47 } 48 49 // -- Get Pre-requisite things 50 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size)); 51 52 PetscCall(GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, -1, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL)); 53 PetscCall(CreateBasisFromPlex(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, &basis_grid_aniso)); 54 55 // -- Build RHS operator 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", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data)); 62 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE)); 63 64 PetscCall(OperatorApplyContextCreate(user->dm, grid_aniso_proj->dm, ceed, op_rhs_assemble, CEED_VECTOR_NONE, NULL, NULL, NULL, &l2_rhs_ctx)); 65 66 // -- Build Mass Operator 67 PetscCall(CreateMassQFunction(ceed, grid_aniso_proj->num_comp, q_data_size, &qf_mass)); 68 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 69 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE)); 70 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data)); 71 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE)); 72 73 { // -- Setup KSP for L^2 projection 74 Mat mat_mass; 75 PetscCall(OperatorApplyContextCreate(grid_aniso_proj->dm, grid_aniso_proj->dm, ceed, op_mass, NULL, NULL, NULL, NULL, &mass_matop_ctx)); 76 PetscCall(CreateMatShell_Ceed(mass_matop_ctx, &mat_mass)); 77 78 PetscCall(KSPCreate(comm, &ksp)); 79 PetscCall(KSPSetOptionsPrefix(ksp, "grid_anisotropy_tensor_projection_")); 80 { 81 PC pc; 82 PetscCall(KSPGetPC(ksp, &pc)); 83 PetscCall(PCSetType(pc, PCJACOBI)); 84 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 85 PetscCall(KSPSetType(ksp, KSPCG)); 86 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 87 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 88 } 89 PetscCall(KSPSetOperators(ksp, mat_mass, mat_mass)); 90 PetscCall(KSPSetFromOptions(ksp)); 91 } 92 93 { // -- Project anisotropy data and store in CeedVector 94 Vec Grid_Anisotropy, grid_anisotropy_loc; 95 96 // Get L^2 Projection RHS 97 PetscCall(DMGetGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 98 99 PetscCall(ApplyCeedOperatorLocalToGlobal(NULL, Grid_Anisotropy, l2_rhs_ctx)); 100 101 // Solve projection problem 102 PetscCall(KSPSolve(ksp, Grid_Anisotropy, Grid_Anisotropy)); 103 104 // Copy anisotropy tensor data to CeedVector 105 PetscCall(DMGetLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 106 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, grid_aniso_vector, NULL)); 107 PetscCall(DMGlobalToLocal(grid_aniso_proj->dm, Grid_Anisotropy, INSERT_VALUES, grid_anisotropy_loc)); 108 PetscCall(VecCopyP2C(grid_anisotropy_loc, *grid_aniso_vector)); 109 PetscCall(DMRestoreLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc)); 110 PetscCall(DMRestoreGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy)); 111 } 112 113 // -- Cleanup 114 PetscCall(NodalProjectionDataDestroy(grid_aniso_proj)); 115 PetscCall(OperatorApplyContextDestroy(l2_rhs_ctx)); 116 PetscCall(OperatorApplyContextDestroy(mass_matop_ctx)); 117 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_assemble)); 118 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 119 PetscCallCeed(ceed, CeedBasisDestroy(&basis_grid_aniso)); 120 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_assemble)); 121 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 122 PetscCall(KSPDestroy(&ksp)); 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso, 127 CeedVector *aniso_colloc_ceed, PetscInt *num_comp_aniso) { 128 CeedInt q_data_size, num_nodes; 129 CeedQFunction qf_colloc; 130 CeedOperator op_colloc; 131 132 PetscFunctionBeginUser; 133 *num_comp_aniso = 7; 134 PetscCallCeed(ceed, CeedBasisGetNumNodes(ceed_data->basis_q, &num_nodes)); 135 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size)); 136 PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_nodes, *num_comp_aniso, NULL, NULL, elem_restr_grid_aniso)); 137 138 // -- Build collocation operator 139 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc)); 140 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_colloc, "qdata", q_data_size, CEED_EVAL_NONE)); 141 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_colloc, "v", *num_comp_aniso, CEED_EVAL_NONE)); 142 143 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_colloc, NULL, NULL, &op_colloc)); 144 PetscCallCeed(ceed, CeedOperatorSetField(op_colloc, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data)); 145 PetscCallCeed(ceed, CeedOperatorSetField(op_colloc, "v", *elem_restr_grid_aniso, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE)); 146 147 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, aniso_colloc_ceed, NULL)); 148 149 PetscCallCeed(ceed, CeedOperatorApply(op_colloc, CEED_VECTOR_NONE, *aniso_colloc_ceed, CEED_REQUEST_IMMEDIATE)); 150 151 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_colloc)); 152 PetscCallCeed(ceed, CeedOperatorDestroy(&op_colloc)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155