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 */ 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 36 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, label_value, height, 1, elem_restr_inv_multiplicity)); 37 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_inv_multiplicity, inv_multiplicity, NULL)); 38 39 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, dm_field, &elem_restr)); 40 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr, &num_comp)); 41 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr, &multiplicity, NULL)); 42 43 if (get_global_multiplicity) { 44 // In order to get global multiplicity, we need to run through DMLocalToGlobal -> DMGlobalToLocal. 45 PetscCall(DMGetLocalVector(dm, &Multiplicity_loc)); 46 PetscCall(DMGetGlobalVector(dm, &Multiplicity)); 47 48 PetscCall(VecPetscToCeed(Multiplicity_loc, &m_mem_type, multiplicity)); 49 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr, multiplicity)); 50 PetscCall(VecCeedToPetsc(multiplicity, m_mem_type, Multiplicity_loc)); 51 PetscCall(VecZeroEntries(Multiplicity)); 52 PetscCall(DMLocalToGlobal(dm, Multiplicity_loc, ADD_VALUES, Multiplicity)); 53 PetscCall(DMGlobalToLocal(dm, Multiplicity, INSERT_VALUES, Multiplicity_loc)); 54 PetscCall(VecPetscToCeed(Multiplicity_loc, &m_mem_type, multiplicity)); 55 } else { 56 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr, multiplicity)); 57 } 58 59 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, InverseMultiplicity, InverseMultiplicity_loc, &qf_multiplicity)); 60 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_multiplicity, "multiplicity", num_comp, CEED_EVAL_NONE)); 61 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_multiplicity, "inverse multiplicity", 1, CEED_EVAL_NONE)); 62 63 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_multiplicity, NULL, NULL, &op_multiplicity)); 64 PetscCallCeed(ceed, CeedOperatorSetName(op_multiplicity, "InverseMultiplicity")); 65 PetscCallCeed(ceed, CeedOperatorSetField(op_multiplicity, "multiplicity", elem_restr, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 66 PetscCallCeed(ceed, 67 CeedOperatorSetField(op_multiplicity, "inverse multiplicity", *elem_restr_inv_multiplicity, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 68 69 PetscCallCeed(ceed, CeedOperatorApply(op_multiplicity, multiplicity, *inv_multiplicity, CEED_REQUEST_IMMEDIATE)); 70 71 if (get_global_multiplicity) { 72 PetscCall(VecCeedToPetsc(multiplicity, m_mem_type, Multiplicity_loc)); 73 PetscCall(DMRestoreLocalVector(dm, &Multiplicity_loc)); 74 PetscCall(DMRestoreGlobalVector(dm, &Multiplicity)); 75 } 76 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity)); 77 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr)); 78 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_multiplicity)); 79 PetscCallCeed(ceed, CeedOperatorDestroy(&op_multiplicity)); 80 PetscFunctionReturn(PETSC_SUCCESS); 81 } 82