PetscDiffSolver and adaptivity

[salazardetroya]

I modified example fem_system4 to use PetscDiffSolver instead and I ran into two issues:

First, the solver is destroyed after each solve(), so when we call equation_systems.reinit(), petsc cannot find the KSP. This can be fixed with a simple PetscDiffSolver::init() in reinit()

Second, there is a lack a proper convergence after an AMR. The problem is linear, but it is not being solved in one iteration even with direct solvers. Before the AMR, it is solved with just one iteration. I am going to try to investigate this and send a patch.

This is the modified fem_system.C https://paste.ofcode.org/TNEuaCWUTPbUqyuTEEQygG and this is the modified heat_system.C https://paste.ofcode.org/33SgJeFhZaZHKdB5jJZnK49

Running with -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps -ksp_converged_reason -ksp_final_residual -log_summary -snes_monitor -ksp_view -snes_monitor_residual. I see the ksp final residual is pretty much zero, but not the snes residual. (EDIT: Just realized the ksp final residual does not have to coincide with the snes residual...)

I also run the 2x1 simple case (no AMR) with -snes_type test -snes_test_display shows an error in the jacobian, but it looks like it is just a failure to apply the BC for the first two degrees of freedom in the residual.

Finite difference Jacobian (constant state 1.0) Mat Object: 1 MPI processes type: seqaij row 0: (0, 0.) row 1: (1, 0.) row 2: (0, 0.333333) (1, 0.166667) (2, -1.33333) (3, 0.333333) (4, 0.166667) (5, 0.333333) row 3: (0, 0.166667) (1, 0.333333) (2, 0.333333) (3, -1.33333) (4, 0.333333) (5, 0.166667) row 4: (2, 0.166667) (3, 0.333333) (4, -0.666667) (5, 0.166667) row 5: (2, 0.333333) (3, 0.166667) (4, 0.166667) (5, -0.666667) Hand-coded Jacobian (constant state 1.0) Mat Object:() 1 MPI processes type: seqaij row 0: (0, 1.) (1, 0.) (2, 0.) (3, 0.) row 1: (0, 0.) (1, 1.) (2, 0.) (3, 0.) row 2: (0, 0.) (1, 0.) (2, -1.33333) (3, 0.333333) (4, 0.166667) (5, 0.333333) row 3: (0, 0.) (1, 0.) (2, 0.333333) (3, -1.33333) (4, 0.333333) (5, 0.166667) row 4: (2, 0.166667) (3, 0.333333) (4, -0.666667) (5, 0.166667) row 5: (2, 0.333333) (3, 0.166667) (4, 0.166667) (5, -0.666667) Hand-coded minus finite-difference Jacobian (constant state 1.0) Mat Object: 1 MPI processes type: seqaij row 0: (0, 1.) (1, 0.) (2, 0.) (3, 0.) row 1: (0, 0.) (1, 1.) (2, 0.) (3, 0.) row 2: (0, -0.333333) (1, -0.166667) (2, -9.39529e-09) (3, 1.43835e-08) (4, 3.18019e-09) (5, 6.36038e-09) row 3: (0, -0.166667) (1, -0.333333) (2, -9.68584e-09) (3, -1.74184e-08) (4, -1.66273e-09) (5, -4.84292e-09) row 4: (2, 3.18019e-09) (3, -1.66273e-09) (4, -4.69764e-09) (5, 3.18019e-09) row 5: (2, -1.66273e-09) (3, -4.84292e-09) (4, 3.18019e-09) (5, -4.69764e-09) Norm of matrix ratio 0.559444, difference 1.50923 (constant state 1.0)

[salazardetroya]

I solved this issue. It has to do with how the hanging node constraints are enforced (https://github.com/libMesh/libmesh/blob/master/src/solvers/petsc_diff_solver.C#L180) the I will submit a PR soon.

[pbauman]

I'd like to keep this issue open until you've posted the PR you mentioned, just so we don't forget in case your efforts towards a PR (which are very welcome!) get interrupted. Thanks!

[ghost]

Ok, this issue is more complicated that it seemed (this is salazardetroya, decided to open a new account).

The line 180 in petsc_diff_solver.C is useless because we are calling System::update() afterwards, returning the local vectors to the state before DofMap::enforce_constraints_exactly(). Therefore, it should be called before enforcing the constraints. However, System::update() is called again inside of ImplicitSystem::assembly(), this should be removed if we want DofMap::enforce_constraints_exactly() to have any effect, but probably this would break other parts of the code that call ImplicitSystem::assembly() so I don't know how to go about this. Any suggestions? Maybe a flag for update in the argument list of System::assembly()? Or maybe move DofMap::enforce_constraints_exactly() inside of ImplicitSystem::assembly() and call it after System::update().