ceed-preconditioning.c (fddff3480f915e14d0e1949ec881b24caa5f6f46) ceed-preconditioning.c (edb2538e3dd6743c029967fc4e89c6fcafedb8c2)
1// Copyright (c) 2017-2022, 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 <ceed-impl.h>

--- 744 unchanged lines hidden (view full) ---

753
754 @param[in] op_fine Fine grid operator
755 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
756 @param[in] rstr_coarse Coarse grid restriction
757 @param[in] basis_coarse Coarse grid active vector basis
758 @param[in] basis_c_to_f Basis for coarse to fine interpolation, or NULL if not creating prolongation/restriction operators
759 @param[out] op_coarse Coarse grid operator
760 @param[out] op_prolong Coarse to fine operator, or NULL
1// Copyright (c) 2017-2022, 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 <ceed-impl.h>

--- 744 unchanged lines hidden (view full) ---

753
754 @param[in] op_fine Fine grid operator
755 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
756 @param[in] rstr_coarse Coarse grid restriction
757 @param[in] basis_coarse Coarse grid active vector basis
758 @param[in] basis_c_to_f Basis for coarse to fine interpolation, or NULL if not creating prolongation/restriction operators
759 @param[out] op_coarse Coarse grid operator
760 @param[out] op_prolong Coarse to fine operator, or NULL
761 @param[out] op_restrict Fine to coarse operator, or NULL
761 @param[out] op_rstrict Fine to coarse operator, or NULL
762
763 @return An error code: 0 - success, otherwise - failure
764
765 @ref Developer
766**/
767static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
762
763 @return An error code: 0 - success, otherwise - failure
764
765 @ref Developer
766**/
767static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
768 CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) {
768 CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_rstrict) {
769 bool is_composite;
770 Ceed ceed;
771 CeedInt num_comp;
772 CeedVector mult_vec = NULL;
773 CeedElemRestriction rstr_p_mult_fine = NULL, rstr_fine = NULL;
774
775 CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
776

--- 21 unchanged lines hidden (view full) ---

798 CeedCall(CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name, op_fine->output_fields[i]->elem_rstr,
799 op_fine->output_fields[i]->basis, op_fine->output_fields[i]->vec));
800 }
801 }
802 // -- Clone QFunctionAssemblyData
803 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(op_fine->qf_assembled, &(*op_coarse)->qf_assembled));
804
805 // Multiplicity vector
769 bool is_composite;
770 Ceed ceed;
771 CeedInt num_comp;
772 CeedVector mult_vec = NULL;
773 CeedElemRestriction rstr_p_mult_fine = NULL, rstr_fine = NULL;
774
775 CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
776

--- 21 unchanged lines hidden (view full) ---

798 CeedCall(CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name, op_fine->output_fields[i]->elem_rstr,
799 op_fine->output_fields[i]->basis, op_fine->output_fields[i]->vec));
800 }
801 }
802 // -- Clone QFunctionAssemblyData
803 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(op_fine->qf_assembled, &(*op_coarse)->qf_assembled));
804
805 // Multiplicity vector
806 if (op_restrict || op_prolong) {
806 if (op_rstrict || op_prolong) {
807 CeedVector mult_e_vec;
808 CeedRestrictionType rstr_type;
809
810 CeedCall(CeedElemRestrictionGetType(rstr_fine, &rstr_type));
811 CeedCheck(rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_UNSUPPORTED,
812 "Element restrictions created with CeedElemRestrictionCreateCurlOriented are not supported");
813 CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector");
814 CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_fine));

--- 6 unchanged lines hidden (view full) ---

821 CeedCall(CeedVectorReciprocal(mult_vec));
822 }
823
824 // Clone name
825 bool has_name = op_fine->name;
826 size_t name_len = op_fine->name ? strlen(op_fine->name) : 0;
827 CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name));
828
807 CeedVector mult_e_vec;
808 CeedRestrictionType rstr_type;
809
810 CeedCall(CeedElemRestrictionGetType(rstr_fine, &rstr_type));
811 CeedCheck(rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_UNSUPPORTED,
812 "Element restrictions created with CeedElemRestrictionCreateCurlOriented are not supported");
813 CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector");
814 CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_fine));

--- 6 unchanged lines hidden (view full) ---

821 CeedCall(CeedVectorReciprocal(mult_vec));
822 }
823
824 // Clone name
825 bool has_name = op_fine->name;
826 size_t name_len = op_fine->name ? strlen(op_fine->name) : 0;
827 CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name));
828
829 // Check that coarse to fine basis is provided if prolong/restrict operators are requested
830 CeedCheck(basis_c_to_f || (!op_restrict && !op_prolong), ceed, CEED_ERROR_INCOMPATIBLE,
829 // Check that coarse to fine basis is provided if prolong/rstrict operators are requested
830 CeedCheck(basis_c_to_f || (!op_rstrict && !op_prolong), ceed, CEED_ERROR_INCOMPATIBLE,
831 "Prolongation or restriction operator creation requires coarse-to-fine basis");
832
833 // Restriction/Prolongation Operators
834 CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp));
835
836 // Restriction
831 "Prolongation or restriction operator creation requires coarse-to-fine basis");
832
833 // Restriction/Prolongation Operators
834 CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp));
835
836 // Restriction
837 if (op_restrict) {
837 if (op_rstrict) {
838 CeedInt *num_comp_r_data;
839 CeedQFunctionContext ctx_r;
838 CeedInt *num_comp_r_data;
839 CeedQFunctionContext ctx_r;
840 CeedQFunction qf_restrict;
840 CeedQFunction qf_rstrict;
841
841
842 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_restrict));
842 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_rstrict));
843 CeedCall(CeedCalloc(1, &num_comp_r_data));
844 num_comp_r_data[0] = num_comp;
845 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r));
846 CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data));
843 CeedCall(CeedCalloc(1, &num_comp_r_data));
844 num_comp_r_data[0] = num_comp;
845 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r));
846 CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data));
847 CeedCall(CeedQFunctionSetContext(qf_restrict, ctx_r));
847 CeedCall(CeedQFunctionSetContext(qf_rstrict, ctx_r));
848 CeedCall(CeedQFunctionContextDestroy(&ctx_r));
848 CeedCall(CeedQFunctionContextDestroy(&ctx_r));
849 CeedCall(CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE));
850 CeedCall(CeedQFunctionAddInput(qf_restrict, "scale", num_comp, CEED_EVAL_NONE));
851 CeedCall(CeedQFunctionAddOutput(qf_restrict, "output", num_comp, CEED_EVAL_INTERP));
852 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_restrict, num_comp));
849 CeedCall(CeedQFunctionAddInput(qf_rstrict, "input", num_comp, CEED_EVAL_NONE));
850 CeedCall(CeedQFunctionAddInput(qf_rstrict, "scale", num_comp, CEED_EVAL_NONE));
851 CeedCall(CeedQFunctionAddOutput(qf_rstrict, "output", num_comp, CEED_EVAL_INTERP));
852 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_rstrict, num_comp));
853
853
854 CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict));
855 CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
856 CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec));
857 CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE));
854 CeedCall(CeedOperatorCreate(ceed, qf_rstrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_rstrict));
855 CeedCall(CeedOperatorSetField(*op_rstrict, "input", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
856 CeedCall(CeedOperatorSetField(*op_rstrict, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec));
857 CeedCall(CeedOperatorSetField(*op_rstrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE));
858
859 // Set name
860 char *restriction_name;
861
862 CeedCall(CeedCalloc(17 + name_len, &restriction_name));
863 sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : "");
858
859 // Set name
860 char *restriction_name;
861
862 CeedCall(CeedCalloc(17 + name_len, &restriction_name));
863 sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : "");
864 CeedCall(CeedOperatorSetName(*op_restrict, restriction_name));
864 CeedCall(CeedOperatorSetName(*op_rstrict, restriction_name));
865 CeedCall(CeedFree(&restriction_name));
866
867 // Check
865 CeedCall(CeedFree(&restriction_name));
866
867 // Check
868 CeedCall(CeedOperatorCheckReady(*op_restrict));
868 CeedCall(CeedOperatorCheckReady(*op_rstrict));
869
870 // Cleanup
869
870 // Cleanup
871 CeedCall(CeedQFunctionDestroy(&qf_restrict));
871 CeedCall(CeedQFunctionDestroy(&qf_rstrict));
872 }
873
874 // Prolongation
875 if (op_prolong) {
876 CeedInt *num_comp_p_data;
877 CeedQFunctionContext ctx_p;
878 CeedQFunction qf_prolong;
879

--- 295 unchanged lines hidden (view full) ---

1175**/
1176int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) {
1177 CeedInt num_active_bases = 0, num_input_fields, *num_eval_modes_in = NULL, *num_eval_modes_out = NULL, offset = 0, num_output_fields;
1178 CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL;
1179 CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL;
1180 CeedQFunctionField *qf_fields;
1181 CeedQFunction qf;
1182 CeedOperatorField *op_fields;
872 }
873
874 // Prolongation
875 if (op_prolong) {
876 CeedInt *num_comp_p_data;
877 CeedQFunctionContext ctx_p;
878 CeedQFunction qf_prolong;
879

--- 295 unchanged lines hidden (view full) ---

1175**/
1176int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) {
1177 CeedInt num_active_bases = 0, num_input_fields, *num_eval_modes_in = NULL, *num_eval_modes_out = NULL, offset = 0, num_output_fields;
1178 CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL;
1179 CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL;
1180 CeedQFunctionField *qf_fields;
1181 CeedQFunction qf;
1182 CeedOperatorField *op_fields;
1183 bool is_composite;
1184
1183
1185 CeedCall(CeedOperatorIsComposite(op, &is_composite));
1186 CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Can only create CeedOperator assembly data for non-composite operators.");
1187
1188 // Allocate
1189 CeedCall(CeedCalloc(1, data));
1190 (*data)->ceed = ceed;
1191 CeedCall(CeedReference(ceed));
1192
1193 // Build OperatorAssembly data
1194 CeedCall(CeedOperatorGetQFunction(op, &qf));
1195 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL));

--- 582 unchanged lines hidden (view full) ---

1778 CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, false, assembled));
1779 } else {
1780 CeedCall(CeedSingleOperatorAssembleAddDiagonal_Core(op, request, false, assembled));
1781 }
1782 return CEED_ERROR_SUCCESS;
1783}
1784
1785/**
1184 // Allocate
1185 CeedCall(CeedCalloc(1, data));
1186 (*data)->ceed = ceed;
1187 CeedCall(CeedReference(ceed));
1188
1189 // Build OperatorAssembly data
1190 CeedCall(CeedOperatorGetQFunction(op, &qf));
1191 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL));

--- 582 unchanged lines hidden (view full) ---

1774 CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, false, assembled));
1775 } else {
1776 CeedCall(CeedSingleOperatorAssembleAddDiagonal_Core(op, request, false, assembled));
1777 }
1778 return CEED_ERROR_SUCCESS;
1779}
1780
1781/**
1786 @brief Fully assemble the point-block diagonal pattern of a linear operator.
1787
1788 Expected to be used in conjunction with CeedOperatorLinearAssemblePointBlockDiagonal().
1789
1790 The assembly routines use coordinate format, with `num_entries` tuples of the form (i, j, value) which indicate that value should be added to the
1791matrix in entry (i, j).
1792 Note that the (i, j) pairs are unique.
1793 This function returns the number of entries and their (i, j) locations, while CeedOperatorLinearAssemblePointBlockDiagonal() provides the values in
1794the same ordering.
1795
1796 This will generally be slow unless your operator is low-order.
1797
1798 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable.
1799
1800 @param[in] op CeedOperator to assemble
1801 @param[out] num_entries Number of entries in coordinate nonzero pattern
1802 @param[out] rows Row number for each entry
1803 @param[out] cols Column number for each entry
1804
1805 @ref User
1806**/
1807int CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) {
1808 Ceed ceed;
1809 bool is_composite;
1810 CeedInt num_active_components, num_sub_operators;
1811 CeedOperator *sub_operators;
1812
1813 CeedCall(CeedOperatorGetCeed(op, &ceed));
1814 CeedCall(CeedOperatorIsComposite(op, &is_composite));
1815
1816 CeedSize input_size = 0, output_size = 0;
1817 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
1818 CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square");
1819
1820 if (is_composite) {
1821 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub_operators));
1822 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
1823 } else {
1824 sub_operators = &op;
1825 num_sub_operators = 1;
1826 }
1827
1828 { // Verify operator can be assembled correctly
1829 CeedInt num_active_elem_rstrs, comp_stride;
1830 CeedOperatorAssemblyData data;
1831 CeedElemRestriction *active_elem_rstrs;
1832
1833 // Get initial values to check against
1834 CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[0], &data));
1835 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs));
1836 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[0], &comp_stride));
1837 CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[0], &num_active_components));
1838
1839 for (CeedInt k = 0; k < num_sub_operators; k++) {
1840 CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[k], &data));
1841
1842 // Verify that all active element restrictions have same component stride and number of components
1843 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs));
1844 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[0], &comp_stride));
1845 for (CeedInt i = 0; i < num_active_elem_rstrs; i++) {
1846 CeedInt comp_stride_sub;
1847 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[i], &comp_stride_sub));
1848 CeedCheck(comp_stride == comp_stride_sub, ceed, CEED_ERROR_DIMENSION,
1849 "Active element restrictions must have the same component stride: %d vs %d", comp_stride, comp_stride_sub);
1850
1851 CeedInt num_active_components_sub;
1852 CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[i], &num_active_components_sub));
1853 CeedCheck(num_active_components == num_active_components_sub, ceed, CEED_ERROR_INCOMPATIBLE,
1854 "All suboperators must have the same number of output components");
1855 }
1856 }
1857 }
1858
1859 *num_entries = input_size * num_active_components;
1860 CeedCall(CeedCalloc(*num_entries, rows));
1861 CeedCall(CeedCalloc(*num_entries, cols));
1862
1863 for (CeedInt o = 0; o < num_sub_operators; o++) {
1864 CeedElemRestriction active_elem_rstr, pb_active_elem_rstr;
1865 CeedInt comp_stride, num_elem, elem_size;
1866 const CeedInt *offsets, *pb_offsets;
1867
1868 CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[o], &active_elem_rstr));
1869 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstr, &comp_stride));
1870 CeedCall(CeedElemRestrictionGetNumElements(active_elem_rstr, &num_elem));
1871 CeedCall(CeedElemRestrictionGetElementSize(active_elem_rstr, &elem_size));
1872 CeedCall(CeedElemRestrictionGetOffsets(active_elem_rstr, CEED_MEM_HOST, &offsets));
1873
1874 CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstr, &pb_active_elem_rstr));
1875 CeedCall(CeedElemRestrictionGetOffsets(pb_active_elem_rstr, CEED_MEM_HOST, &pb_offsets));
1876
1877 for (CeedSize i = 0; i < num_elem * elem_size; i++) {
1878 for (CeedInt c_out = 0; c_out < num_active_components; c_out++) {
1879 for (CeedInt c_in = 0; c_in < num_active_components; c_in++) {
1880 (*rows)[pb_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_out * comp_stride;
1881 (*cols)[pb_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_in * comp_stride;
1882 }
1883 }
1884 }
1885
1886 CeedCall(CeedElemRestrictionRestoreOffsets(active_elem_rstr, &offsets));
1887 CeedCall(CeedElemRestrictionRestoreOffsets(pb_active_elem_rstr, &pb_offsets));
1888 CeedCall(CeedElemRestrictionDestroy(&pb_active_elem_rstr));
1889 }
1890 return CEED_ERROR_SUCCESS;
1891}
1892
1893/**
1894 @brief Assemble the point block diagonal of a square linear CeedOperator
1895
1896 This overwrites a CeedVector with the point block diagonal of a linear CeedOperator.
1897
1898 Note: Currently only non-composite CeedOperators with a single field and composite CeedOperators with single field sub-operators are supported.
1899
1900 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable.
1901

--- 335 unchanged lines hidden (view full) ---

2237 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable.
2238
2239 @param[in] op_fine Fine grid operator
2240 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
2241 @param[in] rstr_coarse Coarse grid restriction
2242 @param[in] basis_coarse Coarse grid active vector basis
2243 @param[out] op_coarse Coarse grid operator
2244 @param[out] op_prolong Coarse to fine operator, or NULL
1782 @brief Assemble the point block diagonal of a square linear CeedOperator
1783
1784 This overwrites a CeedVector with the point block diagonal of a linear CeedOperator.
1785
1786 Note: Currently only non-composite CeedOperators with a single field and composite CeedOperators with single field sub-operators are supported.
1787
1788 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable.
1789

--- 335 unchanged lines hidden (view full) ---

2125 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable.
2126
2127 @param[in] op_fine Fine grid operator
2128 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
2129 @param[in] rstr_coarse Coarse grid restriction
2130 @param[in] basis_coarse Coarse grid active vector basis
2131 @param[out] op_coarse Coarse grid operator
2132 @param[out] op_prolong Coarse to fine operator, or NULL
2245 @param[out] op_restrict Fine to coarse operator, or NULL
2133 @param[out] op_rstrict Fine to coarse operator, or NULL
2246
2247 @return An error code: 0 - success, otherwise - failure
2248
2249 @ref User
2250**/
2251int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
2134
2135 @return An error code: 0 - success, otherwise - failure
2136
2137 @ref User
2138**/
2139int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
2252 CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) {
2140 CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_rstrict) {
2253 CeedBasis basis_c_to_f = NULL;
2254
2255 CeedCall(CeedOperatorCheckReady(op_fine));
2256
2257 // Build prolongation matrix, if required
2141 CeedBasis basis_c_to_f = NULL;
2142
2143 CeedCall(CeedOperatorCheckReady(op_fine));
2144
2145 // Build prolongation matrix, if required
2258 if (op_prolong || op_restrict) {
2146 if (op_prolong || op_rstrict) {
2259 CeedBasis basis_fine;
2260
2261 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
2262 CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f));
2263 }
2264
2265 // Core code
2147 CeedBasis basis_fine;
2148
2149 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
2150 CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f));
2151 }
2152
2153 // Core code
2266 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict));
2154 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_rstrict));
2267 return CEED_ERROR_SUCCESS;
2268}
2269
2270/**
2271 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator with a tensor basis for the active basis
2272
2273 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable.
2274
2275 @param[in] op_fine Fine grid operator
2276 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
2277 @param[in] rstr_coarse Coarse grid restriction
2278 @param[in] basis_coarse Coarse grid active vector basis
2279 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or NULL if not creating prolongation/restriction operators
2280 @param[out] op_coarse Coarse grid operator
2281 @param[out] op_prolong Coarse to fine operator, or NULL
2155 return CEED_ERROR_SUCCESS;
2156}
2157
2158/**
2159 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator with a tensor basis for the active basis
2160
2161 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable.
2162
2163 @param[in] op_fine Fine grid operator
2164 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
2165 @param[in] rstr_coarse Coarse grid restriction
2166 @param[in] basis_coarse Coarse grid active vector basis
2167 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or NULL if not creating prolongation/restriction operators
2168 @param[out] op_coarse Coarse grid operator
2169 @param[out] op_prolong Coarse to fine operator, or NULL
2282 @param[out] op_restrict Fine to coarse operator, or NULL
2170 @param[out] op_rstrict Fine to coarse operator, or NULL
2283
2284 @return An error code: 0 - success, otherwise - failure
2285
2286 @ref User
2287**/
2288int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
2289 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong,
2171
2172 @return An error code: 0 - success, otherwise - failure
2173
2174 @ref User
2175**/
2176int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
2177 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong,
2290 CeedOperator *op_restrict) {
2178 CeedOperator *op_rstrict) {
2291 Ceed ceed;
2292 CeedInt Q_f, Q_c;
2293 CeedBasis basis_fine, basis_c_to_f = NULL;
2294
2295 CeedCall(CeedOperatorCheckReady(op_fine));
2296 CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
2297
2298 // Check for compatible quadrature spaces
2299 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
2300 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f));
2301 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c));
2302 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces");
2303
2304 // Create coarse to fine basis, if required
2179 Ceed ceed;
2180 CeedInt Q_f, Q_c;
2181 CeedBasis basis_fine, basis_c_to_f = NULL;
2182
2183 CeedCall(CeedOperatorCheckReady(op_fine));
2184 CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
2185
2186 // Check for compatible quadrature spaces
2187 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
2188 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f));
2189 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c));
2190 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces");
2191
2192 // Create coarse to fine basis, if required
2305 if (op_prolong || op_restrict) {
2193 if (op_prolong || op_rstrict) {
2306 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c;
2307 CeedScalar *q_ref, *q_weight, *grad;
2308
2309 // Check if interpolation matrix is provided
2310 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE,
2311 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix");
2312 CeedCall(CeedBasisGetDimension(basis_fine, &dim));
2313 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp));

--- 5 unchanged lines hidden (view full) ---

2319 CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad));
2320 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_c, P_1d_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f));
2321 CeedCall(CeedFree(&q_ref));
2322 CeedCall(CeedFree(&q_weight));
2323 CeedCall(CeedFree(&grad));
2324 }
2325
2326 // Core code
2194 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c;
2195 CeedScalar *q_ref, *q_weight, *grad;
2196
2197 // Check if interpolation matrix is provided
2198 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE,
2199 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix");
2200 CeedCall(CeedBasisGetDimension(basis_fine, &dim));
2201 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp));

--- 5 unchanged lines hidden (view full) ---

2207 CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad));
2208 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_c, P_1d_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f));
2209 CeedCall(CeedFree(&q_ref));
2210 CeedCall(CeedFree(&q_weight));
2211 CeedCall(CeedFree(&grad));
2212 }
2213
2214 // Core code
2327 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict));
2215 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_rstrict));
2328 return CEED_ERROR_SUCCESS;
2329}
2330
2331/**
2332 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator with a non-tensor basis for the active vector
2333
2334 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable.
2335
2336 @param[in] op_fine Fine grid operator
2337 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
2338 @param[in] rstr_coarse Coarse grid restriction
2339 @param[in] basis_coarse Coarse grid active vector basis
2340 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or NULL if not creating prolongation/restriction operators
2341 @param[out] op_coarse Coarse grid operator
2342 @param[out] op_prolong Coarse to fine operator, or NULL
2216 return CEED_ERROR_SUCCESS;
2217}
2218
2219/**
2220 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator with a non-tensor basis for the active vector
2221
2222 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable.
2223
2224 @param[in] op_fine Fine grid operator
2225 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or NULL if not creating prolongation/restriction operators
2226 @param[in] rstr_coarse Coarse grid restriction
2227 @param[in] basis_coarse Coarse grid active vector basis
2228 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or NULL if not creating prolongation/restriction operators
2229 @param[out] op_coarse Coarse grid operator
2230 @param[out] op_prolong Coarse to fine operator, or NULL
2343 @param[out] op_restrict Fine to coarse operator, or NULL
2231 @param[out] op_rstrict Fine to coarse operator, or NULL
2344
2345 @return An error code: 0 - success, otherwise - failure
2346
2347 @ref User
2348**/
2349int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
2232
2233 @return An error code: 0 - success, otherwise - failure
2234
2235 @ref User
2236**/
2237int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
2350 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong,
2351 CeedOperator *op_restrict) {
2238 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_rstrict) {
2352 Ceed ceed;
2353 CeedInt Q_f, Q_c;
2354 CeedBasis basis_fine, basis_c_to_f = NULL;
2355
2356 CeedCall(CeedOperatorCheckReady(op_fine));
2357 CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
2358
2359 // Check for compatible quadrature spaces
2360 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
2361 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f));
2362 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c));
2363 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces");
2364
2365 // Coarse to fine basis
2239 Ceed ceed;
2240 CeedInt Q_f, Q_c;
2241 CeedBasis basis_fine, basis_c_to_f = NULL;
2242
2243 CeedCall(CeedOperatorCheckReady(op_fine));
2244 CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
2245
2246 // Check for compatible quadrature spaces
2247 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
2248 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f));
2249 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c));
2250 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces");
2251
2252 // Coarse to fine basis
2366 if (op_prolong || op_restrict) {
2253 if (op_prolong || op_rstrict) {
2367 CeedInt dim, num_comp, num_nodes_c, num_nodes_f;
2368 CeedScalar *q_ref, *q_weight, *grad;
2369 CeedElemTopology topo;
2370
2371 // Check if interpolation matrix is provided
2372 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE,
2373 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix");
2374 CeedCall(CeedBasisGetTopology(basis_fine, &topo));

--- 6 unchanged lines hidden (view full) ---

2381 CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad));
2382 CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_c, num_nodes_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f));
2383 CeedCall(CeedFree(&q_ref));
2384 CeedCall(CeedFree(&q_weight));
2385 CeedCall(CeedFree(&grad));
2386 }
2387
2388 // Core code
2254 CeedInt dim, num_comp, num_nodes_c, num_nodes_f;
2255 CeedScalar *q_ref, *q_weight, *grad;
2256 CeedElemTopology topo;
2257
2258 // Check if interpolation matrix is provided
2259 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE,
2260 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix");
2261 CeedCall(CeedBasisGetTopology(basis_fine, &topo));

--- 6 unchanged lines hidden (view full) ---

2268 CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad));
2269 CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_c, num_nodes_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f));
2270 CeedCall(CeedFree(&q_ref));
2271 CeedCall(CeedFree(&q_weight));
2272 CeedCall(CeedFree(&grad));
2273 }
2274
2275 // Core code
2389 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict));
2276 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_rstrict));
2390 return CEED_ERROR_SUCCESS;
2391}
2392
2393/**
2394 @brief Build a FDM based approximate inverse for each element for a CeedOperator
2395
2396 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization Method based approximate inverse.
2397 This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, \f$M = V^T V, K = V^T S V\f$.

--- 234 unchanged lines hidden --- 		2277 return CEED_ERROR_SUCCESS;
2278}
2279
2280/**
2281 @brief Build a FDM based approximate inverse for each element for a CeedOperator
2282
2283 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization Method based approximate inverse.
2284 This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, \f$M = V^T V, K = V^T S V\f$.

--- 234 unchanged lines hidden ---