xref: /honee/src/qdata.c (revision 15c15616687dd5d637ad23a3aaa500a53e25a78e) 
1 // Copyright (c) 2017-2024, 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 <navierstokes.h>
9 
10 #include <petscsection.h>
11 #include "../qfunctions/setupgeo.h"
12 #include "../qfunctions/setupgeo2d.h"
13 
14 /**
15  * @brief Get number of components of quadrature data for domain
16  *
17  * @param[in]  dm          DM where quadrature data would be used
18  * @param[out] q_data_size Number of components of quadrature data
19  */
20 PetscErrorCode QDataGetNumComponents(DM dm, CeedInt *q_data_size) {
21   PetscInt num_comp_x, dim;
22 
23   PetscFunctionBeginUser;
24   PetscCall(DMGetDimension(dm, &dim));
25   {  // Get number of coordinate components
26     DM           dm_coord;
27     PetscSection section_coord;
28     PetscInt     field = 0;  // Default field has the coordinates
29     PetscCall(DMGetCoordinateDM(dm, &dm_coord));
30     PetscCall(DMGetLocalSection(dm_coord, &section_coord));
31     PetscCall(PetscSectionGetFieldComponents(section_coord, field, &num_comp_x));
32   }
33   switch (dim) {
34     case 2:
35       switch (num_comp_x) {
36         case 2:
37           *q_data_size = 5;
38           break;
39         case 3:
40           *q_data_size = 7;
41           break;
42         default:
43           SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
44                   "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
45           break;
46       }
47       break;
48     case 3:
49       *q_data_size = 10;
50       break;
51     default:
52       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
53               "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
54       break;
55   }
56   PetscFunctionReturn(PETSC_SUCCESS);
57 }
58 
59 /**
60  * @brief Create quadrature data for domain
61  *
62  * @param[in]  ceed          Ceed object quadrature data will be used with
63  * @param[in]  dm            DM where quadrature data would be used
64  * @param[in]  domain_label  DMLabel that quadrature data would be used one
65  * @param[in]  label_value   Value of label
66  * @param[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
67  * @param[in]  basis_x       CeedBasis of the coordinates
68  * @param[in]  x_coord       CeedVector of the coordinates
69  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
70  * @param[out] q_data        CeedVector of the quadrature data
71  * @param[out] q_data_size   number of components of quadrature data
72  */
73 PetscErrorCode QDataGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
74                         CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
75   CeedQFunction qf_setup = NULL;
76   CeedOperator  op_setup;
77   CeedInt       num_comp_x;
78   PetscInt      dim, height = 0;
79 
80   PetscFunctionBeginUser;
81   PetscCall(QDataGetNumComponents(dm, q_data_size));
82   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
83   PetscCall(DMGetDimension(dm, &dim));
84   switch (dim) {
85     case 2:
86       switch (num_comp_x) {
87         case 2:
88           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2d, Setup2d_loc, &qf_setup));
89           break;
90         case 3:
91           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2D_3Dcoords, Setup2D_3Dcoords_loc, &qf_setup));
92           break;
93       }
94       break;
95     case 3:
96       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup, Setup_loc, &qf_setup));
97       break;
98   }
99 
100   // -- Create QFunction for quadrature data
101   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup, 0));
102   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
103   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT));
104   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "surface qdata", *q_data_size, CEED_EVAL_NONE));
105 
106   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, elem_restr_qd));
107   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_qd, q_data, NULL));
108 
109   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup));
110   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
111   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
112   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "surface qdata", *elem_restr_qd, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
113 
114   PetscCallCeed(ceed, CeedOperatorApply(op_setup, x_coord, *q_data, CEED_REQUEST_IMMEDIATE));
115 
116   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
117   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
118   PetscFunctionReturn(PETSC_SUCCESS);
119 }
120 
121 /**
122  * @brief Get number of components of quadrature data for boundary of domain
123  *
124  * @param[in]  dm          DM where quadrature data would be used
125  * @param[out] q_data_size Number of components of quadrature data
126  */
127 PetscErrorCode QDataBoundaryGetNumComponents(DM dm, CeedInt *q_data_size) {
128   PetscInt dim;
129 
130   PetscFunctionBeginUser;
131   PetscCall(DMGetDimension(dm, &dim));
132   switch (dim) {
133     case 2:
134       *q_data_size = 3;
135       break;
136     case 3:
137       *q_data_size = 10;
138       break;
139     default:
140       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "QDataBoundary not valid for DM of dimension %" PetscInt_FMT, dim);
141       break;
142   }
143   PetscFunctionReturn(PETSC_SUCCESS);
144 }
145 
146 /**
147  * @brief Create quadrature data for boundary of domain
148  *
149  * @param[in]  ceed          Ceed object quadrature data will be used with
150  * @param[in]  dm            DM where quadrature data would be used
151  * @param[in]  domain_label  DMLabel that quadrature data would be used one
152  * @param[in]  label_value   Value of label
153  * @param[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
154  * @param[in]  basis_x       CeedBasis of the coordinates
155  * @param[in]  x_coord       CeedVector of the coordinates
156  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
157  * @param[out] q_data        CeedVector of the quadrature data
158  * @param[out] q_data_size   number of components of quadrature data
159  */
160 PetscErrorCode QDataBoundaryGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
161                                 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
162   CeedQFunction qf_setup_sur = NULL;
163   CeedOperator  op_setup_sur;
164   CeedInt       num_comp_x;
165   PetscInt      dim, height = 1;
166 
167   PetscFunctionBeginUser;
168   PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size));
169   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
170   PetscCall(DMGetDimension(dm, &dim));
171   switch (dim) {
172     case 2:
173       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary2d, SetupBoundary2d_loc, &qf_setup_sur));
174       break;
175     case 3:
176       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary, SetupBoundary_loc, &qf_setup_sur));
177       break;
178   }
179 
180   // -- Create QFunction for quadrature data
181   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup_sur, 0));
182   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
183   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
184   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "surface qdata", *q_data_size, CEED_EVAL_NONE));
185 
186   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, elem_restr_qd));
187   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_qd, q_data, NULL));
188 
189   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, NULL, NULL, &op_setup_sur));
190   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
191   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
192   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "surface qdata", *elem_restr_qd, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
193 
194   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, *q_data, CEED_REQUEST_IMMEDIATE));
195 
196   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
197   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
198   PetscFunctionReturn(PETSC_SUCCESS);
199 }
200