xref: /honee/src/grid_anisotropy_tensor.c (revision 00359db47665a79ecb0241f6ccbf886b649022df) 
1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 
4 #include "../qfunctions/grid_anisotropy_tensor.h"
5 
6 #include <petscdmplex.h>
7 
8 #include <navierstokes.h>
9 
GridAnisotropyTensorProjectionSetupApply(Ceed ceed,Honee honee,CeedElemRestriction * elem_restr_grid_aniso,CeedVector * grid_aniso_vector)10 PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, Honee honee, CeedElemRestriction *elem_restr_grid_aniso,
11                                                         CeedVector *grid_aniso_vector) {
12   NodalProjectionData grid_aniso_proj;
13   CeedBasis           basis_grid_aniso;
14   CeedVector          q_data;
15   CeedElemRestriction elem_restr_qd;
16   CeedInt             q_data_size;
17   MPI_Comm            comm   = PetscObjectComm((PetscObject)honee->dm);
18   PetscInt            height = 0, dm_field = 0;
19 
20   PetscFunctionBeginUser;
21   PetscCall(PetscNew(&grid_aniso_proj));
22 
23   {  // -- Create DM for Anisotropic tensor L^2 projection
24     PetscSection section;
25 
26     PetscCall(DMClone(honee->dm, &grid_aniso_proj->dm));
27     PetscCall(DMSetMatrixPreallocateSkip(grid_aniso_proj->dm, PETSC_TRUE));
28     PetscCall(PetscObjectSetName((PetscObject)grid_aniso_proj->dm, "Grid Anisotropy Tensor Projection"));
29 
30     // -- Setup DM
31     grid_aniso_proj->num_comp = 7;
32     PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, honee->app_ctx->degree, 1, honee->app_ctx->q_extra, 1, &grid_aniso_proj->num_comp,
33                                  grid_aniso_proj->dm));
34 
35     PetscCall(DMGetLocalSection(grid_aniso_proj->dm, &section));
36     PetscCall(PetscSectionSetFieldName(section, 0, ""));
37     PetscCall(PetscSectionSetComponentName(section, 0, 0, "KMGridAnisotropyTensorXX"));
38     PetscCall(PetscSectionSetComponentName(section, 0, 1, "KMGridAnisotropyTensorYY"));
39     PetscCall(PetscSectionSetComponentName(section, 0, 2, "KMGridAnisotropyTensorZZ"));
40     PetscCall(PetscSectionSetComponentName(section, 0, 3, "KMGridAnisotropyTensorYZ"));
41     PetscCall(PetscSectionSetComponentName(section, 0, 4, "KMGridAnisotropyTensorXZ"));
42     PetscCall(PetscSectionSetComponentName(section, 0, 5, "KMGridAnisotropyTensorXY"));
43     PetscCall(PetscSectionSetComponentName(section, 0, 6, "GridAnisotropyTensorFrobNorm"));
44   }
45 
46   // -- Get Pre-requisite things
47   PetscCall(DMPlexCeedElemRestrictionCreate(ceed, grid_aniso_proj->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field,
48                                             elem_restr_grid_aniso));
49   PetscCall(DMPlexCeedBasisCreate(ceed, grid_aniso_proj->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, &basis_grid_aniso));
50   PetscCall(QDataGet(ceed, grid_aniso_proj->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size));
51 
52   {  // -- Build RHS operator
53     CeedOperator  op_rhs_assemble;
54     CeedQFunction qf_rhs_assemble;
55 
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", elem_restr_qd, CEED_BASIS_NONE, q_data));
62     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE));
63 
64     PetscCall(OperatorApplyContextCreate(honee->dm, grid_aniso_proj->dm, ceed, op_rhs_assemble, CEED_VECTOR_NONE, NULL, NULL, NULL,
65                                          &grid_aniso_proj->l2_rhs_ctx));
66 
67     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_assemble));
68     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_assemble));
69   }
70 
71   {  // Setup KSP for L^2 projection
72     CeedOperator  op_mass;
73     CeedQFunction qf_mass;
74     Mat           mat_mass;
75 
76     PetscCall(HoneeMassQFunctionCreate(ceed, grid_aniso_proj->num_comp, q_data_size, &qf_mass));
77     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
78     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE));
79     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
80     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE));
81 
82     PetscCall(MatCreateCeed(grid_aniso_proj->dm, grid_aniso_proj->dm, op_mass, NULL, &mat_mass));
83 
84     PetscCall(KSPCreate(comm, &grid_aniso_proj->ksp));
85     PetscCall(KSPSetOptionsPrefix(grid_aniso_proj->ksp, "grid_anisotropy_tensor_projection_"));
86     {
87       PC pc;
88       PetscCall(KSPGetPC(grid_aniso_proj->ksp, &pc));
89       PetscCall(PCSetType(pc, PCJACOBI));
90       PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL));
91       PetscCall(KSPSetType(grid_aniso_proj->ksp, KSPCG));
92       PetscCall(KSPSetNormType(grid_aniso_proj->ksp, KSP_NORM_NATURAL));
93       PetscCall(KSPSetTolerances(grid_aniso_proj->ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
94     }
95     PetscCall(KSPSetFromOptions_WithMatCeed(grid_aniso_proj->ksp, mat_mass));
96 
97     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
98     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
99     PetscCall(MatDestroy(&mat_mass));
100   }
101 
102   {  // -- Project anisotropy data and store in CeedVector
103     Vec Grid_Anisotropy, grid_anisotropy_loc;
104 
105     PetscCall(DMGetGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy));
106     PetscCall(ApplyCeedOperatorLocalToGlobal(NULL, Grid_Anisotropy, grid_aniso_proj->l2_rhs_ctx));
107     PetscCall(KSPSolve(grid_aniso_proj->ksp, Grid_Anisotropy, Grid_Anisotropy));
108 
109     // Copy anisotropy tensor data to CeedVector
110     PetscCall(DMGetLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc));
111     PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, grid_aniso_vector, NULL));
112     PetscCall(DMGlobalToLocal(grid_aniso_proj->dm, Grid_Anisotropy, INSERT_VALUES, grid_anisotropy_loc));
113     PetscCall(VecCopyPetscToCeed(grid_anisotropy_loc, *grid_aniso_vector));
114     PetscCall(DMRestoreLocalVector(grid_aniso_proj->dm, &grid_anisotropy_loc));
115     PetscCall(DMRestoreGlobalVector(grid_aniso_proj->dm, &Grid_Anisotropy));
116   }
117 
118   PetscCall(NodalProjectionDataDestroy(&grid_aniso_proj));
119   PetscCallCeed(ceed, CeedBasisDestroy(&basis_grid_aniso));
120   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
121   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
122   PetscFunctionReturn(PETSC_SUCCESS);
123 }
124 
GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed,Honee honee,CeedElemRestriction * elem_restr_grid_aniso,CeedVector * aniso_colloc_ceed,PetscInt * num_comp_aniso)125 PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, Honee honee, CeedElemRestriction *elem_restr_grid_aniso,
126                                                              CeedVector *aniso_colloc_ceed, PetscInt *num_comp_aniso) {
127   CeedInt             q_data_size, num_nodes;
128   CeedQFunction       qf_colloc;
129   CeedOperator        op_colloc;
130   CeedVector          q_data;
131   CeedElemRestriction elem_restr_qd;
132   PetscInt            height = 0;
133 
134   PetscFunctionBeginUser;
135   *num_comp_aniso = 7;
136   {
137     CeedBasis basis_q;
138     PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q));
139     PetscCallCeed(ceed, CeedBasisGetNumNodes(basis_q, &num_nodes));
140     PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
141   }
142   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, *num_comp_aniso,
143                                                  elem_restr_grid_aniso));
144   PetscCall(QDataGet(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size));
145 
146   // -- Build collocation operator
147   PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc));
148   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_colloc, "qdata", q_data_size, CEED_EVAL_NONE));
149   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_colloc, "v", *num_comp_aniso, CEED_EVAL_NONE));
150 
151   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_colloc, NULL, NULL, &op_colloc));
152   PetscCallCeed(ceed, CeedOperatorSetField(op_colloc, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
153   PetscCallCeed(ceed, CeedOperatorSetField(op_colloc, "v", *elem_restr_grid_aniso, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
154 
155   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, aniso_colloc_ceed, NULL));
156 
157   PetscCallCeed(ceed, CeedOperatorApply(op_colloc, CEED_VECTOR_NONE, *aniso_colloc_ceed, CEED_REQUEST_IMMEDIATE));
158 
159   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
160   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
161   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_colloc));
162   PetscCallCeed(ceed, CeedOperatorDestroy(&op_colloc));
163   PetscFunctionReturn(PETSC_SUCCESS);
164 }
165