Forms the complete direct FD Jacobian, checks every row, and compares MFFD actions.
665{
668 SNES snes = NULL;
669 Mat J = NULL;
670 Vec x = NULL, xp = NULL, f0 = NULL, fp = NULL, column = NULL, square = NULL;
671 Vec row_norm_sq = NULL, v[2] = {NULL, NULL}, dense_v[2] = {NULL, NULL};
672 Vec mffd_v = NULL, error_vec = NULL;
673 PetscInt n_global, lo, hi;
674 PetscReal min_row_sq = 0.0, error = 0.0, reference = 0.0;
675 const PetscReal epsilon = 1.0e-7;
677
678 PetscFunctionBeginUser;
680 PetscCall(VecDuplicate(user->
Ucont, &x)); PetscCall(VecSet(x, 0.0));
681 PetscCall(VecDuplicate(x, &xp)); PetscCall(VecDuplicate(x, &f0));
682 PetscCall(VecDuplicate(x, &fp)); PetscCall(VecDuplicate(x, &column));
683 PetscCall(VecDuplicate(x, &square)); PetscCall(VecDuplicate(x, &row_norm_sq));
684 PetscCall(VecDuplicate(x, &v[0])); PetscCall(VecDuplicate(x, &v[1]));
685 PetscCall(VecDuplicate(x, &dense_v[0])); PetscCall(VecDuplicate(x, &dense_v[1]));
686 PetscCall(VecDuplicate(x, &mffd_v)); PetscCall(VecDuplicate(x, &error_vec));
687 PetscCall(VecZeroEntries(row_norm_sq)); PetscCall(VecZeroEntries(dense_v[0]));
688 PetscCall(VecZeroEntries(dense_v[1]));
689 PetscCall(VecGetSize(x, &n_global));
690 PetscCall(VecGetOwnershipRange(x, &lo, &hi));
691 for (PetscInt which = 0; which < 2; ++which) {
692 PetscScalar *a = NULL;
693 PetscCall(VecGetArray(v[which], &a));
694 for (PetscInt local = 0; local < hi - lo; ++local) {
695 PetscInt global = lo + local;
696 a[local] = which == 0 ? (PetscScalar)(1.0 + 0.05 * (global % 9))
697 : (PetscScalar)(((global % 2) ? -1.0 : 1.0) * (0.5 + 0.03 * (global % 7)));
698 }
699 PetscCall(VecRestoreArray(v[which], &a));
700 }
701
705 for (PetscInt col = 0; col < n_global; ++col) {
706 PetscScalar coeff[2] = {
707 (PetscScalar)(1.0 + 0.05 * (col % 9)),
708 (PetscScalar)(((col % 2) ? -1.0 : 1.0) * (0.5 + 0.03 * (col % 7)))
709 };
710 PetscCall(VecCopy(x, xp));
711 if (col >= lo && col < hi) PetscCall(VecSetValue(xp, col, epsilon, ADD_VALUES));
712 PetscCall(VecAssemblyBegin(xp)); PetscCall(VecAssemblyEnd(xp));
714 PetscCall(VecWAXPY(column, -1.0, f0, fp));
715 PetscCall(VecScale(column, 1.0 / epsilon));
716 PetscCall(VecPointwiseMult(square, column, column));
717 PetscCall(VecAXPY(row_norm_sq, 1.0, square));
718 PetscCall(VecAXPY(dense_v[0], coeff[0], column));
719 PetscCall(VecAXPY(dense_v[1], coeff[1], column));
720 }
721 PetscCall(VecMin(row_norm_sq, NULL, &min_row_sq));
722 PetscCheck(min_row_sq > 0.5, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
723 "Complete Newton Jacobian contains an unexplained zero/weak row (min squared norm=%g).",
724 (double)min_row_sq);
725 PetscCall(
CheckStoredDerivative(user, x, 2, 2, 2, 0, 2, 2, 2, 0,
726 10.0, 1.0e-5, "interior BDF1 temporal diagonal"));
727
728 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
729 PetscCall(SNESSetDM(snes, user->
fda));
731 PetscCall(MatCreateSNESMF(snes, &J));
732 PetscCall(MatMFFDSetBase(J, x, f0));
733 for (PetscInt which = 0; which < 2; ++which) {
734 PetscCall(MatMult(J, v[which], mffd_v));
735 PetscCall(VecWAXPY(error_vec, -1.0, dense_v[which], mffd_v));
736 PetscCall(VecNorm(error_vec, NORM_2, &error));
737 PetscCall(VecNorm(dense_v[which], NORM_2, &reference));
738 PetscCheck(error <= 2.0e-5 * PetscMax(1.0, reference), PETSC_COMM_WORLD, PETSC_ERR_PLIB,
739 "Independent dense FD Jv differs from PETSc MFFD action %d: error=%g reference=%g.",
740 which, (double)error, (double)reference);
741 }
742
743 for (PetscInt which = 0; which < 2; ++which) {
744 const PetscReal steps[3] = {1.0e-4, 1.0e-6, 1.0e-8};
745 PetscReal best = PETSC_MAX_REAL;
746 for (PetscInt s = 0; s < 3; ++s) {
747 PetscCall(VecWAXPY(xp, steps[s], v[which], x));
749 PetscCall(VecWAXPY(column, -1.0, f0, fp));
750 PetscCall(VecScale(column, 1.0 / steps[s]));
751 PetscCall(VecAXPY(column, -1.0, dense_v[which]));
752 PetscCall(VecNorm(column, NORM_2, &error));
753 best = PetscMin(best, error);
754 }
755 PetscCall(VecNorm(dense_v[which], NORM_2, &reference));
756 PetscCheck(best <= 2.0e-5 * PetscMax(1.0, reference), PETSC_COMM_WORLD, PETSC_ERR_PLIB,
757 "Direct directional differences show no accuracy plateau for vector %d (best=%g).",
758 which, (double)best);
759 }
760
761 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes));
762 PetscCall(VecDestroy(&error_vec)); PetscCall(VecDestroy(&mffd_v));
763 PetscCall(VecDestroy(&dense_v[1])); PetscCall(VecDestroy(&dense_v[0]));
764 PetscCall(VecDestroy(&v[1])); PetscCall(VecDestroy(&v[0]));
765 PetscCall(VecDestroy(&row_norm_sq)); PetscCall(VecDestroy(&square));
766 PetscCall(VecDestroy(&column)); PetscCall(VecDestroy(&fp)); PetscCall(VecDestroy(&f0));
767 PetscCall(VecDestroy(&xp)); PetscCall(VecDestroy(&x));
770 PetscFunctionReturn(PETSC_SUCCESS);
771}
PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
Cleans up and destroys all boundary system resources.
static PetscErrorCode BuildMinimalWallOperatorFixture(SimCtx **simCtx, UserCtx **user, PetscBool x_periodic)
Builds a compact all-wall operator fixture through real boundary handlers.
PetscErrorCode PicurvDestroyMinimalContexts(SimCtx **simCtx_ptr, UserCtx **user_ptr)
Destroys minimal SimCtx/UserCtx fixtures and all owned PETSc objects.
