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