1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 #include "../qfunctions/inverse_multiplicity.h"
5 #include <navierstokes.h>
6
7 /**
8 * @brief Get the inverse of the node multiplicity
9 *
10 * Local multiplicity concerns only the number of elements associated with a node in the current rank's partition.
11 * Global multiplicity is the multiplicity for the global mesh, including periodicity.
12 *
13 * @param[in] ceed `Ceed` object for the output `CeedVector` and `CeedElemRestriction`
14 * @param[in] dm `DM` for the grid
15 * @param[in] domain_label `DMLabel` for `DMPlex` domain
16 * @param[in] label_value Stratum value
17 * @param[in] height Height of `DMPlex` topology
18 * @param[in] dm_field Index of `DMPlex` field
19 * @param[in] get_global_multiplicity Whether the multiplicity should be global or local
20 * @param[out] elem_restr_inv_multiplicity `CeedElemRestriction` needed to use the multiplicity vector
21 * @param[out] inv_multiplicity `CeedVector` containing the inverse of the multiplicity
22 */
GetInverseMultiplicity(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt dm_field,PetscBool get_global_multiplicity,CeedElemRestriction * elem_restr_inv_multiplicity,CeedVector * inv_multiplicity)23 PetscErrorCode GetInverseMultiplicity(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field,
24 PetscBool get_global_multiplicity, CeedElemRestriction *elem_restr_inv_multiplicity,
25 CeedVector *inv_multiplicity) {
26 Vec Multiplicity, Multiplicity_loc;
27 PetscMemType m_mem_type;
28 CeedVector multiplicity;
29 CeedQFunction qf_multiplicity;
30 CeedOperator op_multiplicity;
31 CeedInt num_comp;
32 CeedElemRestriction elem_restr;
33
34 PetscFunctionBeginUser;
35 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, label_value, height, 1, elem_restr_inv_multiplicity));
36 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_inv_multiplicity, inv_multiplicity, NULL));
37
38 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, dm_field, &elem_restr));
39 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr, &num_comp));
40 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr, &multiplicity, NULL));
41
42 if (get_global_multiplicity) {
43 // In order to get global multiplicity, we need to run through DMLocalToGlobal -> DMGlobalToLocal.
44 PetscCall(DMGetLocalVector(dm, &Multiplicity_loc));
45 PetscCall(DMGetGlobalVector(dm, &Multiplicity));
46
47 PetscCall(VecPetscToCeed(Multiplicity_loc, &m_mem_type, multiplicity));
48 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr, multiplicity));
49 PetscCall(VecCeedToPetsc(multiplicity, m_mem_type, Multiplicity_loc));
50 PetscCall(VecZeroEntries(Multiplicity));
51 PetscCall(DMLocalToGlobal(dm, Multiplicity_loc, ADD_VALUES, Multiplicity));
52 PetscCall(DMGlobalToLocal(dm, Multiplicity, INSERT_VALUES, Multiplicity_loc));
53 PetscCall(VecPetscToCeed(Multiplicity_loc, &m_mem_type, multiplicity));
54 } else {
55 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr, multiplicity));
56 }
57
58 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, InverseMultiplicity, InverseMultiplicity_loc, &qf_multiplicity));
59 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_multiplicity, "multiplicity", num_comp, CEED_EVAL_NONE));
60 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_multiplicity, "inverse multiplicity", 1, CEED_EVAL_NONE));
61
62 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_multiplicity, NULL, NULL, &op_multiplicity));
63 PetscCallCeed(ceed, CeedOperatorSetName(op_multiplicity, "InverseMultiplicity"));
64 PetscCallCeed(ceed, CeedOperatorSetField(op_multiplicity, "multiplicity", elem_restr, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
65 PetscCallCeed(ceed,
66 CeedOperatorSetField(op_multiplicity, "inverse multiplicity", *elem_restr_inv_multiplicity, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
67
68 PetscCallCeed(ceed, CeedOperatorApply(op_multiplicity, multiplicity, *inv_multiplicity, CEED_REQUEST_IMMEDIATE));
69
70 if (get_global_multiplicity) {
71 PetscCall(VecCeedToPetsc(multiplicity, m_mem_type, Multiplicity_loc));
72 PetscCall(DMRestoreLocalVector(dm, &Multiplicity_loc));
73 PetscCall(DMRestoreGlobalVector(dm, &Multiplicity));
74 }
75 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity));
76 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr));
77 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_multiplicity));
78 PetscCallCeed(ceed, CeedOperatorDestroy(&op_multiplicity));
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81