Verifies the point-block model and common preconditioning-engine wiring.
1660{
1663 char tmpdir[PETSC_MAX_PATH_LEN] = "";
1668 0,
1669 0
1670 };
1672 Mat preconditioning_matrix = NULL;
1673 Vec x = NULL, f = NULL;
1674 PetscInt block_size = 0, velocity_dof = 0;
1675 PetscReal matrix_norm = 0.0, reassembled_norm = 0.0, difference_norm = 0.0;
1676 PetscScalar reference[9], values[9], mutant[9], legacy[9];
1677 const PetscInt target_i = 1, target_j = 2, target_k = 3;
1678 MatStencil target[3] = {
1679 {.i = 1, .j = 2, .k = 3, .c = 0},
1680 {.i = 1, .j = 2, .k = 3, .c = 1},
1681 {.i = 1, .j = 2, .k = 3, .c = 2}
1682 };
1683 MatStencil conditioned = {.i = 0, .j = 2, .k = 3, .c = 0};
1684 MatStencil homogeneous = {.i = 0, .j = 2, .k = 3, .c = 1};
1685 Mat saved_matrix = NULL, mffd = NULL;
1686 SNES snes = NULL;
1687 Vec direction = NULL, product = NULL, px = NULL;
1688 PetscScalar px_sum = 0.0;
1689 PetscReal px_norm = 0.0;
1690 KSP ksp = NULL;
1691 PC pc = NULL;
1692 Vec pc_rhs = NULL, pc_solution = NULL;
1693 MatInfo initial_allocation_info, repeated_allocation_info;
1694
1695 PetscFunctionBeginUser;
1697 PetscCall(VecDuplicate(user->
Ucont, &user->
Rhs));
1698 PetscCall(VecDuplicate(user->
Ucont, &x));
1699 PetscCall(VecDuplicate(user->
Ucont, &f));
1700 PetscCall(VecCopy(user->
Ucont, x));
1706 "engine must select the frozen-momentum model callbacks"));
1708 "point-block engine must own its separate matrix"));
1710 "point-block matrix must not alias the Jacobian operator"));
1713 user, preconditioning_matrix, PETSC_FALSE));
1714 PetscCall(MatGetInfo(preconditioning_matrix, MAT_GLOBAL_SUM,
1715 &initial_allocation_info));
1716 PetscCall(MatNorm(preconditioning_matrix, NORM_FROBENIUS, &matrix_norm));
1718 "model callback and common rows must insert matrix entries"));
1719 PetscCall(MatGetBlockSize(preconditioning_matrix, &block_size));
1721 PetscCall(DMDAGetInfo(user->
fda, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
1722 &velocity_dof, NULL, NULL, NULL, NULL, NULL));
1724
1725
1727 for (PetscInt rr = 0; rr < 3; ++rr)
1728 for (PetscInt cc = 0; cc < 3; ++cc) {
1730 target[rr], target[cc], &values[3 * rr + cc]));
1732 PetscRealPart(values[3 * rr + cc]), 1e-13,
1733 "Cartesian point block entry must match independent oracle"));
1734 }
1735
1736
1737 for (PetscInt cc = 0; cc < 3; ++cc) {
1738 MatStencil col = conditioned; col.c = cc;
1739 PetscScalar conditioned_value = 0.0, homogeneous_value = 0.0;
1741 conditioned, col, &conditioned_value));
1742 col = homogeneous; col.c = cc;
1744 homogeneous, col, &homogeneous_value));
1746 PetscRealPart(conditioned_value), 1e-14, "conditioned row must be exact identity"));
1748 PetscRealPart(homogeneous_value), 1e-14, "homogeneous row must be exact identity"));
1749 }
1750
1751
1752 PetscCall(MatDuplicate(preconditioning_matrix, MAT_COPY_VALUES, &saved_matrix));
1753 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
1754 PetscCall(SNESSetDM(snes, user->
fda));
1756 PetscCall(MatCreateSNESMF(snes, &mffd));
1757 PetscCall(VecDuplicate(x, &direction)); PetscCall(VecDuplicate(x, &product));
1758 PetscCall(VecSet(direction, 0.375));
1760 PetscCall(MatMFFDSetBase(mffd, x, f));
1761 PetscCall(MatMult(mffd, direction, product));
1763 PetscCall(MatAXPY(saved_matrix, -1.0, preconditioning_matrix, SAME_NONZERO_PATTERN));
1764 PetscCall(MatNorm(saved_matrix, NORM_FROBENIUS, &difference_norm));
1766 "assembly must be unchanged after residual and MFFD products"));
1767 PetscCall(VecDestroy(&product)); PetscCall(VecDestroy(&direction));
1768 PetscCall(MatDestroy(&mffd)); PetscCall(SNESDestroy(&snes)); PetscCall(MatDestroy(&saved_matrix));
1769
1770
1772 PetscCall(VecCopy(user->
Ucont, x));
1776 for (PetscInt rr = 0; rr < 3; ++rr)
1777 for (PetscInt cc = 0; cc < 3; ++cc) {
1779 target[rr], target[cc], &values[3 * rr + cc]));
1781 PetscRealPart(values[3 * rr + cc]), 1e-13,
1782 "nonuniform point block entry must match independent oracle"));
1783 }
1784 PetscCall(
PicurvAssertBool((PetscBool)(PetscAbsScalar(reference[1] - reference[3]) > 1e-6 &&
1785 PetscAbsScalar(reference[2] - reference[6]) > 1e-6 &&
1786 PetscAbsScalar(reference[5] - reference[7]) > 1e-6),
1787 "oracle must preserve nonsymmetric component ordering"));
1789 mutant_flag <<= 1) {
1790 PetscBool differs = PETSC_FALSE;
1792 mutant_flag, mutant));
1793 for (PetscInt n = 0; n < 9; ++n)
1794 if (PetscAbsScalar(reference[n] - mutant[n]) > 1e-8) differs = PETSC_TRUE;
1795 PetscCall(
PicurvAssertBool(differs,
"independent oracle must reject audited mutant"));
1796 }
1799 PetscCall(
PicurvAssertBool((PetscBool)(PetscAbsScalar(reference[8] - mutant[8]) > 1e-8),
1800 "A[5] must contribute to the zeta diagonal"));
1803 for (PetscInt n = 0; n < 9; ++n)
1805 1e-13, "modern block must be the negative legacy block"));
1806
1807
1809 {
1810 PetscBool differs = PETSC_FALSE;
1811 for (PetscInt n = 0; n < 9; ++n)
1812 if (PetscAbsScalar(reference[n] - mutant[n]) > 1e-8) differs = PETSC_TRUE;
1813 PetscCall(
PicurvAssertBool(differs,
"permuted MatStencil coordinates must be detectable"));
1814 }
1815 {
1816 MatStencil neighbor = target[0];
1817 PetscScalar neighbor_value = 0.0;
1818 neighbor.i++;
1820 target[0], neighbor, &neighbor_value));
1822 "point block must not insert unintended neighbor entries"));
1823 }
1824
1825
1829 for (PetscInt n = 0; n < 9; ++n) {
1831 target[n / 3], target[n % 3], &values[n]));
1833 1e-13, "BDF2 point block must match independent oracle"));
1834 if (n == 0 || n == 4 || n == 8)
1836 PetscRealPart(mutant[n] - reference[n]), 1e-10, "BDF2 diagonal increment"));
1837 else
1839 1e-13, "BDF order must not change off-diagonal entries"));
1840 }
1841
1842 PetscCall(VecDuplicate(x, &px));
1843 PetscCall(MatMult(preconditioning_matrix, x, px));
1844 PetscCall(VecSum(px, &px_sum)); PetscCall(VecNorm(px, NORM_2, &px_norm));
1846 "P*x global sum must be decomposition independent"));
1848 "P*x norm must be decomposition independent"));
1849 PetscCall(PetscPrintf(PETSC_COMM_WORLD,
1850 "POINT_BLOCK_MPI_SIGNATURE sum=%.16e norm2=%.16e\n",
1851 (double)PetscRealPart(px_sum), (double)px_norm));
1852 PetscCall(VecDestroy(&px));
1853
1854 {
1855 PetscBool assembled = PETSC_FALSE;
1856 PetscCall(MatAssembled(preconditioning_matrix, &assembled));
1857 PetscCall(
PicurvAssertBool(assembled,
"engine must perform final matrix assembly"));
1858 }
1859 PetscCall(MatNorm(preconditioning_matrix, NORM_FROBENIUS, &matrix_norm));
1860 PetscCall(MatShift(preconditioning_matrix, 7.0));
1862 PetscCall(MatGetInfo(preconditioning_matrix, MAT_GLOBAL_SUM,
1863 &repeated_allocation_info));
1865 (PetscReal)repeated_allocation_info.nz_allocated, 0.0,
1866 "repeated assembly must retain exact allocated storage"));
1868 (PetscReal)repeated_allocation_info.mallocs, 0.0,
1869 "repeated assembly must not add insertion reallocations"));
1870 PetscCall(MatNorm(preconditioning_matrix, NORM_FROBENIUS, &reassembled_norm));
1872 "repeated engine assembly must clear old entries"));
1873 {
1874 PetscScalar value = 0.0;
1876 conditioned, conditioned, &value));
1878 "repeated engine assembly must clear old entries"));
1879 }
1880 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
1881 PetscCall(KSPSetOperators(ksp, preconditioning_matrix, preconditioning_matrix));
1882 PetscCall(KSPGetPC(ksp, &pc));
1885 PetscCall(KSPSetUp(ksp));
1886 PetscCall(DMCreateGlobalVector(user->
fda, &pc_rhs));
1887 PetscCall(DMCreateGlobalVector(user->
fda, &pc_solution));
1888 PetscCall(VecSet(pc_rhs, 1.0));
1889 PetscCall(PCApply(pc, pc_rhs, pc_solution));
1890 PetscCall(VecDestroy(&pc_solution)); PetscCall(VecDestroy(&pc_rhs));
1891 PetscCall(KSPDestroy(&ksp));
1894 "engine destroy must clear its owned matrix"));
1895 PetscCall(VecDestroy(&f)); PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&user->
Rhs));
1897 PetscFunctionReturn(PETSC_SUCCESS);
1898}
static const MomentumPreconditionerModelOps frozen_momentum_point_block_ops
PetscBool aliases_jacobian_operator
const MomentumPreconditionerModelOps * model_ops
static PetscErrorCode MomentumPreconditionerEngine_ConfigurePetscPC(MomentumPreconditionerEngine *engine, PC pc)
Applies the validated model/structure-to-PETSc-PC mapping.
static PetscErrorCode MomentumPreconditionerEngine_ValidatePetscPC(MomentumPreconditionerEngine *engine, PC pc)
Rejects raw options that select an unvalidated PETSc PC backend.
static PetscErrorCode MomentumPreconditionerEngine_Create(UserCtx *user, Mat jacobian_operator, const MomentumPreconditionerDescription *requested, MomentumPreconditionerEngine *engine)
Validates a model/structure and creates or aliases its matrix.
static PetscErrorCode MomentumPreconditionerEngine_Assemble(MomentumPreconditionerEngine *engine, UserCtx *user, Vec current_solution)
Runs model insertion, common row handling, and final assembly.
Mat preconditioning_matrix
PetscBool owns_preconditioning_matrix
static PetscErrorCode MomentumPreconditionerEngine_Destroy(MomentumPreconditionerEngine *engine)
Destroys only a separately owned preconditioning matrix.
static PetscErrorCode PreconditionerMatrixStencilEntry(UserCtx *user, Mat preconditioning_matrix, MatStencil row, MatStencil col, PetscScalar *value)
Reads one DMDA-stencil matrix entry through collective basis vectors.
static PetscErrorCode CollectiveLegacyPointBlockOracle(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscInt flags, PetscScalar block[9])
Evaluates the independent oracle on the unique owner and broadcasts it.
static PetscErrorCode SeedPointBlockOracleFields(UserCtx *user)
Seeds nonuniform, index-distinguishing coefficient fields.
static PetscErrorCode AssertExactPointBlockMatrixAllocation(UserCtx *user, Mat matrix, PetscBool require_offrank_periodic)
Verifies the exact AIJ layout and preallocation derived from row classes.
