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 "../include/libceedsetup.h"
9
10 #include <stdio.h>
11
12 #include "../include/libceedsetup.h"
13 #include "../include/petscutils.h"
14
15 // -----------------------------------------------------------------------------
16 // Destroy libCEED operator objects
17 // -----------------------------------------------------------------------------
CeedDataDestroy(CeedInt i,CeedData data)18 PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) {
19 PetscFunctionBeginUser;
20 CeedVectorDestroy(&data->q_data);
21 CeedVectorDestroy(&data->x_ceed);
22 CeedVectorDestroy(&data->y_ceed);
23 CeedBasisDestroy(&data->basis_x);
24 CeedBasisDestroy(&data->basis_u);
25 CeedElemRestrictionDestroy(&data->elem_restr_u);
26 CeedElemRestrictionDestroy(&data->elem_restr_x);
27 CeedElemRestrictionDestroy(&data->elem_restr_u_i);
28 CeedElemRestrictionDestroy(&data->elem_restr_qd_i);
29 CeedQFunctionDestroy(&data->qf_apply);
30 CeedOperatorDestroy(&data->op_apply);
31 if (i > 0) {
32 CeedOperatorDestroy(&data->op_prolong);
33 CeedOperatorDestroy(&data->op_restrict);
34 }
35 PetscCall(PetscFree(data));
36 PetscFunctionReturn(PETSC_SUCCESS);
37 };
38
39 // -----------------------------------------------------------------------------
40 // Set up libCEED for a given degree
41 // -----------------------------------------------------------------------------
SetupLibceedByDegree(DM dm,Ceed ceed,CeedInt degree,CeedInt topo_dim,CeedInt q_extra,PetscInt num_comp_x,PetscInt num_comp_u,PetscInt g_size,PetscInt xl_size,BPData bp_data,CeedData data,PetscBool setup_rhs,PetscBool is_fine_level,CeedVector rhs_ceed,CeedVector * target)42 PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree, CeedInt topo_dim, CeedInt q_extra, PetscInt num_comp_x, PetscInt num_comp_u,
43 PetscInt g_size, PetscInt xl_size, BPData bp_data, CeedData data, PetscBool setup_rhs, PetscBool is_fine_level,
44 CeedVector rhs_ceed, CeedVector *target) {
45 DM dm_coord;
46 Vec coords;
47 const PetscScalar *coord_array;
48 CeedBasis basis_x, basis_u;
49 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i;
50 CeedQFunction qf_setup_geo = NULL, qf_apply = NULL;
51 CeedOperator op_setup_geo, op_apply;
52 CeedVector x_coord, q_data, x_ceed, y_ceed;
53 PetscInt c_start, c_end, num_elem;
54 CeedInt num_qpts, q_data_size = bp_data.q_data_size;
55 CeedScalar R = 1; // radius of the sphere
56 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube
57
58 PetscFunctionBeginUser;
59 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
60
61 // CEED bases
62 PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, bp_data, &basis_x));
63 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u));
64
65 // CEED restrictions
66 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x));
67 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u));
68
69 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end));
70 num_elem = c_end - c_start;
71 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts);
72
73 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, num_comp_u * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_u_i);
74 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, q_data_size * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_qd_i);
75
76 // Element coordinates
77 PetscCall(DMGetCoordinatesLocal(dm, &coords));
78 PetscCall(VecGetArrayRead(coords, &coord_array));
79
80 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL);
81 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coord_array);
82 PetscCall(VecRestoreArrayRead(coords, &coord_array));
83
84 // Create the persistent vectors that will be needed in setup and apply
85 CeedVectorCreate(ceed, q_data_size * num_elem * num_qpts, &q_data);
86 CeedVectorCreate(ceed, xl_size, &x_ceed);
87 CeedVectorCreate(ceed, xl_size, &y_ceed);
88
89 if (is_fine_level) {
90 // Create the QFunction that builds the context data
91 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup_geo);
92 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP);
93 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * topo_dim, CEED_EVAL_GRAD);
94 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT);
95 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE);
96
97 // Create the operator that builds the quadrature data
98 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo);
99 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
100 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
101 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
102 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
103
104 // Setup q_data
105 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE);
106
107 // Set up PDE operator
108 PetscBool is_interp = bp_data.in_mode == CEED_EVAL_INTERP;
109 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? topo_dim : 1;
110 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? topo_dim : 1;
111
112 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply);
113 if (bp_data.in_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) {
114 CeedQFunctionAddInput(qf_apply, "u", num_comp_u, CEED_EVAL_INTERP);
115 CeedQFunctionAddInput(qf_apply, "du", num_comp_u * topo_dim, CEED_EVAL_GRAD);
116 } else {
117 CeedQFunctionAddInput(qf_apply, is_interp ? "u" : "du", num_comp_u * in_scale, bp_data.in_mode);
118 }
119 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE);
120 if (bp_data.out_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) {
121 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u, CEED_EVAL_INTERP);
122 CeedQFunctionAddOutput(qf_apply, "dv", num_comp_u * topo_dim, CEED_EVAL_GRAD);
123 } else {
124 CeedQFunctionAddOutput(qf_apply, is_interp ? "v" : "dv", num_comp_u * out_scale, bp_data.out_mode);
125 }
126
127 // Create the mass or diff operator
128 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply);
129 if (bp_data.in_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) {
130 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
131 CeedOperatorSetField(op_apply, "du", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
132 } else {
133 CeedOperatorSetField(op_apply, is_interp ? "u" : "du", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
134 }
135 CeedOperatorSetField(op_apply, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data);
136 if (bp_data.out_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) {
137 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
138 CeedOperatorSetField(op_apply, "dv", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
139 } else {
140 CeedOperatorSetField(op_apply, is_interp ? "v" : "dv", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
141 }
142
143 // Cleanup
144 CeedQFunctionDestroy(&qf_setup_geo);
145 CeedOperatorDestroy(&op_setup_geo);
146 }
147
148 // Set up RHS if needed
149 if (setup_rhs) {
150 CeedQFunction qf_setup_rhs;
151 CeedOperator op_setup_rhs;
152 CeedVectorCreate(ceed, num_elem * num_qpts * num_comp_u, target);
153 // Create the q-function that sets up the RHS and true solution
154 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs);
155 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP);
156 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE);
157 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, CEED_EVAL_NONE);
158 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP);
159
160 // Create the operator that builds the RHS and true solution
161 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs);
162 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
163 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data);
164 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, CEED_BASIS_NONE, *target);
165 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
166
167 // Set up the libCEED context
168 CeedQFunctionContext ctx_rhs_setup;
169 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup);
170 CeedScalar rhs_setup_data[2] = {R, l};
171 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data);
172 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup);
173 CeedQFunctionContextDestroy(&ctx_rhs_setup);
174
175 // Setup RHS and target
176 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE);
177
178 // Cleanup
179 CeedQFunctionDestroy(&qf_setup_rhs);
180 CeedOperatorDestroy(&op_setup_rhs);
181 }
182 // Cleanup
183 CeedVectorDestroy(&x_coord);
184
185 // Save libCEED data required for level
186 data->basis_x = basis_x;
187 data->basis_u = basis_u;
188 data->elem_restr_x = elem_restr_x;
189 data->elem_restr_u = elem_restr_u;
190 data->elem_restr_u_i = elem_restr_u_i;
191 data->elem_restr_qd_i = elem_restr_qd_i;
192 data->qf_apply = qf_apply;
193 data->op_apply = op_apply;
194 data->q_data = q_data;
195 data->x_ceed = x_ceed;
196 data->y_ceed = y_ceed;
197 data->q_data_size = q_data_size;
198 PetscFunctionReturn(PETSC_SUCCESS);
199 };
200
201 // -----------------------------------------------------------------------------
202 // Setup libCEED level transfer operator objects
203 // -----------------------------------------------------------------------------
CeedLevelTransferSetup(DM dm,Ceed ceed,CeedInt level,CeedInt num_comp_u,CeedData * data,BPData bp_data,Vec fine_mult)204 PetscErrorCode CeedLevelTransferSetup(DM dm, Ceed ceed, CeedInt level, CeedInt num_comp_u, CeedData *data, BPData bp_data, Vec fine_mult) {
205 PetscFunctionBeginUser;
206 // Restriction - Fine to corse
207 CeedOperator op_restrict;
208 // Interpolation - Corse to fine
209 CeedOperator op_prolong;
210 // Coarse grid operator
211 CeedOperator op_apply;
212 // Basis
213 CeedBasis basis_u;
214 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u));
215
216 // ---------------------------------------------------------------------------
217 // Coarse Grid, Prolongation, and Restriction Operators
218 // ---------------------------------------------------------------------------
219 // Create the Operators that compute the prolongation and
220 // restriction between the p-multigrid levels and the coarse grid eval.
221 // ---------------------------------------------------------------------------
222 // Place in libCEED array
223 PetscMemType m_mem_type;
224 PetscCall(VecReadP2C(fine_mult, &m_mem_type, data[level]->x_ceed));
225
226 CeedOperatorMultigridLevelCreate(data[level]->op_apply, data[level]->x_ceed, data[level - 1]->elem_restr_u, basis_u, &op_apply, &op_prolong,
227 &op_restrict);
228
229 // Restore PETSc vector
230 PetscCall(VecReadC2P(data[level]->x_ceed, m_mem_type, fine_mult));
231 PetscCall(VecZeroEntries(fine_mult));
232 // -- Save libCEED data
233 data[level - 1]->op_apply = op_apply;
234 data[level]->op_prolong = op_prolong;
235 data[level]->op_restrict = op_restrict;
236
237 CeedBasisDestroy(&basis_u);
238 PetscFunctionReturn(PETSC_SUCCESS);
239 };
240
SetupErrorOperator(DM dm,Ceed ceed,BPData bp_data,CeedInt topo_dim,PetscInt num_comp_x,PetscInt num_comp_u,CeedOperator * op_error)241 PetscErrorCode SetupErrorOperator(DM dm, Ceed ceed, BPData bp_data, CeedInt topo_dim, PetscInt num_comp_x, PetscInt num_comp_u,
242 CeedOperator *op_error) {
243 DM dm_coord;
244 Vec coords;
245 const PetscScalar *coord_array;
246 CeedBasis basis_x, basis_u;
247 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i;
248 CeedQFunction qf_setup_geo, qf_setup_rhs, qf_error;
249 CeedOperator op_setup_geo, op_setup_rhs;
250 CeedVector x_coord, q_data, target, rhs;
251 PetscInt c_start, c_end, num_elem;
252 CeedInt num_qpts, q_data_size = bp_data.q_data_size;
253 CeedScalar R = 1; // radius of the sphere
254 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube
255
256 PetscFunctionBeginUser;
257 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
258
259 // CEED bases
260 PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, bp_data, &basis_x));
261 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u));
262
263 // CEED restrictions
264 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x));
265 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u));
266
267 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end));
268 num_elem = c_end - c_start;
269 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts);
270
271 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, num_comp_u * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_u_i);
272 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, q_data_size * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_qd_i);
273
274 // Element coordinates
275 PetscCall(DMGetCoordinatesLocal(dm, &coords));
276 PetscCall(VecGetArrayRead(coords, &coord_array));
277
278 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL);
279 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coord_array);
280 PetscCall(VecRestoreArrayRead(coords, &coord_array));
281
282 // Create the persistent vectors that will be needed in setup and apply
283 CeedVectorCreate(ceed, q_data_size * num_elem * num_qpts, &q_data);
284
285 // Create the QFunction that builds the context data
286 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup_geo);
287 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP);
288 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * topo_dim, CEED_EVAL_GRAD);
289 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT);
290 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE);
291
292 // Create the operator that builds the quadrature data
293 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo);
294 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
295 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
296 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
297 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
298
299 // Setup q_data
300 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE);
301
302 // Set up target vector
303 CeedElemRestrictionCreateVector(elem_restr_u, &rhs, NULL);
304 CeedVectorCreate(ceed, num_elem * num_qpts * num_comp_u, &target);
305 // Create the q-function that sets up the RHS and true solution
306 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs);
307 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP);
308 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE);
309 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, CEED_EVAL_NONE);
310 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP);
311
312 // Create the operator that builds the RHS and true solution
313 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs);
314 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
315 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data);
316 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, CEED_BASIS_NONE, target);
317 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
318
319 // Set up the libCEED context
320 CeedQFunctionContext ctx_rhs_setup;
321 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup);
322 CeedScalar rhs_setup_data[2] = {R, l};
323 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data);
324 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup);
325 CeedQFunctionContextDestroy(&ctx_rhs_setup);
326
327 // Setup RHS and target
328 CeedOperatorApply(op_setup_rhs, x_coord, rhs, CEED_REQUEST_IMMEDIATE);
329
330 // Set up error operator
331 // Create the error QFunction
332 CeedQFunctionCreateInterior(ceed, 1, bp_data.error, bp_data.error_loc, &qf_error);
333 CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP);
334 CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE);
335 CeedQFunctionAddInput(qf_error, "qdata", q_data_size, CEED_EVAL_NONE);
336 CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_INTERP);
337
338 // Create the error operator
339 CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_error);
340 CeedOperatorSetField(*op_error, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
341 CeedOperatorSetField(*op_error, "true_soln", elem_restr_u_i, CEED_BASIS_NONE, target);
342 CeedOperatorSetField(*op_error, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data);
343 CeedOperatorSetField(*op_error, "error", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE);
344
345 // Cleanup
346 CeedQFunctionDestroy(&qf_setup_rhs);
347 CeedOperatorDestroy(&op_setup_rhs);
348 CeedQFunctionDestroy(&qf_setup_geo);
349 CeedOperatorDestroy(&op_setup_geo);
350 CeedQFunctionDestroy(&qf_error);
351 CeedVectorDestroy(&x_coord);
352 CeedVectorDestroy(&rhs);
353 CeedVectorDestroy(&target);
354 CeedVectorDestroy(&q_data);
355 CeedElemRestrictionDestroy(&elem_restr_u_i);
356 CeedElemRestrictionDestroy(&elem_restr_qd_i);
357 CeedElemRestrictionDestroy(&elem_restr_x);
358 CeedElemRestrictionDestroy(&elem_restr_u);
359 CeedBasisDestroy(&basis_x);
360 CeedBasisDestroy(&basis_u);
361
362 PetscFunctionReturn(PETSC_SUCCESS);
363 }
364