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