Fix: Update no flux bcs for mesh deformation

[bobmyhill]

This PR is an extension of #6819, and represents a first step towards resolving #6818.

Changes in this PR:

• move the call to compute_no_normal_flux_constraints into compute_current_velocity_boundary_constraints to make sure they are updated during the model run (e.g. if the mesh boundary changes)
• make sure compute_current_velocity_boundary_constraints is called again after mesh deformation if mesh deformation is active to ensure that all boundary conditions are applied correctly to the new (deformed) mesh
• make sure that we do not use the manifold to (re)compute the normals of the deformed mesh.
• adds two new tests.

As described in #6818 this does not yet fix the related GMG issues.

[tjhei]

103 tests failed out of 1279

We will need to take a look at those changes.

[bobmyhill]

@tjhei A few questions:

• Is there a reason that dealii-9.6 has so many more failures than deal.II-9.7 (7, 378) and deal.II-master (7, 676, 680)?
• Of the failures for deal.II-9.7 and deal.II-master, Tests 676 and 680 are related to ghost vector test failure with deal.II master #6822 - I assume we can ignore these in this PR?
• Test 7 works fine on my machine - How can I investigate this failure further?
• Test 378 does appear related to this PR, but implicitly uses GMG. Using AMG allows the test to run (though I haven't yet checked how the output changes). Any thoughts on what I should do with this one?

Test names and errors for deal.II-9.7 and deal.II-master:

• 7 - adiabatic_boundary: An error occurred in line <513> of file </opt/tmp/unpack/deal.II-master/source/matrix_free/dof_info.cc> : Index 183 is not in the half-open range [0,183).
• 378 - geoid_freesurface: Iterative method reported convergence failure in step 1000.
• 676 - nonlinear_channel_flow_velocities_Newton_Stokes: Ghost vector issue (see ghost vector test failure with deal.II master #6822)
• 680 - nonlinear_channel_flow_velocities_Newton_Stokes_no_deviator: Ghost vector issue (see ghost vector test failure with deal.II master #6822)

[tjhei]

It might be worth looking at your and Rene's changes separately.

The deal.II 9.6 tester is the only one that compares test results. The others only execute tests and report crashes. This is because output is slightly different with different deal.II versions.

[tjhei]

I think this is not quite right: If we use mesh deformation we use a flat manifold and we should not use the manifold nornals (as we found out). But without mesh deformation, for example for spherical geometries, we use a different manifold that does provide better normals. So, I would make this bool argument conditional on that.

[bobmyhill]

Hi @tjhei, yes, sorry, I was just trying a few things out to see which change caused the 100-odd tests to fail.

Now that I've pinned down the source of most of the errors, my next step is to use the parameters.mesh_deformation_enabled bool as the argument to the update_no_normal_flux_constraints() calls in

• source/simulator/solver/stokes_matrix_free_global_coarsening.cc
• source/simulator/solver/stokes_matrix_free_local_smoothing.cc
• source/simulator/core.cc

I think this is what you're suggesting?

Any thoughts about what we should do about the failed tests that won't be affected by this change (mesh_deformation_gmg_2d_spherical_shell and four others that use GMG and mesh deformation)?

[tjhei]

The other changes are likely due to Rene's change. If they are small and look sensible, we can update the blessed output.

[bobmyhill]

@tjhei I deleted Rene's change from this PR. The shared commit ("Add test") is just a new test.

The five tests that fail all have significant differences (all involve mesh deformation and the GMG solver):

• 48 - ascii_data_boundary_traction_3d_shell_spherical - Max velocity is significantly different from main (5.01 m/s vs 4.76 m/s). Same for GMG and AMG.
• 122 - boundary_traction_function_cartesian_free_surface - Significant changes to the components of the stress tensor averaged over the domain (>0.1% change)
• 123 - boundary_traction_function_spherical_free_surface - As for 122
• 378 - geoid_freesurface - Now takes a lot more GMG iterations (120-odd vs 90-odd). Other changes seem negligible but maybe the output isn't catching the differences.
• 626 - mesh_deformation_gmg_2d_spherical_shell - The output is significantly different for GMG vs AMG solvers. GMG has a two ?spurious peaks in deformation on either side of the raised elevation

[tjhei]

378 - geoid_freesurface - Now takes a lot more GMG iterations (120-odd vs 90-odd).

... but it might be solving the correct problem now. :-)

626 - mesh_deformation_gmg_2d_spherical_shell - The output is significantly different for GMG vs AMG solvers. GMG has a two ?spurious peaks in deformation on either side of the raised elevation

There are likely other fixes necessary to make GMG work correctly, like you describe in #6818 (comment)

I think this PR is reasonable to merge.

Sorry @tjhei, we'll need to take another look at this - when I removed Rene's first commit and reran the new tests, I forgot to refresh paraview.

I suspect this means I'll have to look through more failing tests, which I'll do tomorrow.

[tiannh7]

Hi @bobmyhill @tjhei,

Thank you very much for this PR! I noticed an interesting behavior during testing and would like to ask for your thoughts.

In the test model ascii_data_boundary_traction_3d_shell_spherical.prm, the top boundary has Free Surface + Traction conditions, while the bottom is Free Slip (no mesh deformation). According to the PR's logic, since mesh deformation is globally enabled, use_manifold_for_normal is set to false, which changes the normal vector computation for the bottom Free Slip boundary as well, ultimately affecting the overall results (e.g., max velocity changes from 4.76 m/s to 5.01 m/s).

My preliminary thought is: Should we check more precisely whether each boundary is actually deformed? For example, only set use_manifold_for_normal = false when the Free Slip boundary itself is deformed. This could avoid introducing unnecessary numerical changes to undeformed boundaries.

What are your opinions on this? Is boundary-specific logic worth considering, or is global consistency more important?

Thanks!

[tjhei]

What are your opinions on this? Is boundary-specific logic worth considering, or is global consistency more important?

Interesting. I didn't think of this situation. Is this an important case worth handling?

We can either determine for each boundary what to do or, maybe the better way, figure out while using the manifold normal doesn't work correctly.

[bobmyhill]

@tiannh7 Thanks for doing some testing!

There are two separate points relating to your questions:

1. This PR is a bug fix for mesh deformation. Is there a reason that you believe the changes to the ascii_data_boundary_traction_3d_shell_spherical output to be due to a regression (a new inaccuracy or bug) rather than the bug fix? I find this test particularly difficult to interrogate.

2. You raise an interesting suggestion for an enhancement to the current PR. It would require a bit of thought - we would either need to loop over all boundaries with tangential velocity boundary conditions to check if they were deforming, or create a new object similar to boundary_velocity_manager.get_tangential_boundary_velocity_indicators() but for only deforming boundaries.

@tjhei: Did you suggest the current logic in this PR because the manifold normals are more accurate / more efficient to calculate?

[tiannh7]

Hi @tjhei and @bobmyhill,

Thank you both for the replies!

Regarding the concern about regression:
My main point is that for static boundaries—like the bottom of a spherical shell—the Manifold Normal provides the exact analytical vector, whereas the Mesh Normal is a discrete approximation. If we disable the Manifold Normal globally, we force static boundaries to use a lower-precision approximation simply because the top boundary is deforming. While this fixes the physical consistency for the deformed surface, it might unintentionally degrade the accuracy on the static boundaries.

I am currently working with several spherical shell models involving free surfaces and traction boundaries. Could I take a little time to run my existing benchmarks against this branch? I want to verify whether the deviation introduced on the static boundaries is negligible or significant. This should help us decide if the current global switch is sufficient or if we typically need finer-grained control per boundary.

When my schedule permits, I’ll report back with my test results as soon as possible!

As a preliminary note: for the ascii_data_boundary_traction_3d_shell_spherical test case, if I apply a patch where use_manifold is set to false only for the deforming boundary (and kept true for static ones), the Maximum velocity is 4.76 m/year.

[tjhei]

@tjhei: Did you suggest the current logic in this PR because the manifold normals are more accurate / more efficient to calculate?

Manifold normals are correct for spherical objects (pointing exactly in radial direction) whereas using the geometry gives slightly different answers, especially on coarse meshes. I think @tiannh7 is right and we need to keep using them for the bottom boundary even if the top is a deforming surface.

And yes, I am suggesting we call compute_no_normal_flux separately for deforming and non-deforming surfaces.

[bobmyhill]

Thanks both. I had a go at making the required change (see most recent commit), but get the error

--------------------------------------------------------
An error occurred in line <891> of file <.../dealii-candi/tmp/unpack/deal.II-v9.7.0/include/deal.II/lac/affine_constraints.templates.h> in function
    void dealii::AffineConstraints<>::close() [number = double]
The violated condition was: 
    line.entries[entry].first != constrained_line.entries[i].first
Additional information: 
    Cycle in constraints detected!

I'm out of time today (I've been doing this in spare 10 minute slots), but any help would be appreciated.

[tiannh7]

Hi Bob,

I believe the following section should be deleted:

VectorTools::compute_no_normal_flux_constraints(dof_handler_v,
                                                      /* first_vector_component= */
                                                      0,
                                                      this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators(),
                                                      constraints_v,
                                                      this->get_mapping(),
                                                      /*use_manifold_for_normal=*/
                                                      !this->get_parameters().mesh_deformation_enabled);

Because it conflicts with sim.compute_current_velocity_boundary_constraints(constraints_v).

[bobmyhill]

Hi @tiannh7
Thanks! I hadn't spotted the duplicated call.

[tiannh7]

@bobmyhill Glad I could help. Looking forward to the next steps.

[bobmyhill]

@tjhei, I think this is ready for another look now.

I reviewed the new tests (which both look fine), and also three of the tests whose output changed in interesting ways:

hill_diffusion_field_boundary_conditions.prm (layer 1)

The Layer 1 composition of the top boundary is determined by the function if(t==0,1,0) if new material is added to the domain. I think the new results look good - there is an overshoot-undershoot artefact, but the transition from p(Layer 1) = 1 to p(Layer 1) = 0 now happens at the correct position (where the boundary is increasing in elevation).

Timestep 0 (both new and old)

End time (new)

End time (old)

boundary_traction_function_cartesian_free_surface.prm (stress xy)

Both solutions have a weird artifact in the middle of the domain. New not obviously worse than old.

new

old

gmg_mesh_deform_topo.prm (velocity and streamlines)

New looks better (streamlines now close). But I did want to note that this test throws up another problem; the shape of the deformed surface at time t is dependent on the size of the maximum time step (converging with a small timestep size). See the two sets of images below.

new (end time, maximum timestep not set)

old (end time, maximum timestep not set)

new (end time, maximum timestep set to 0.00005, global refinement = 5)

old (end time, maximum timestep set to 0.00005, global refinement = 5)

Conclusion

Examples 1 and 3 are clearly improved by this PR, and Example 2 looks much the same (differences are hard to spot over the large unrelated anomalies - Issue #6835 ). The no normal flux constraints now look to be correct for both AMG and GMG, although differences remain in the solutions (which needs further testing, probably outside this PR - Issue #6833). The undershoot problem in Example 1 is also unrelated to this PR (Issue #6834)

[tjhei]

Suggested change

	using the local mesh. This fix focuses on simulations using the
	AMG solver.
	using the mesh.

[bobmyhill]

Thanks, changed!

[tjhei]

Would you mind squashing your commits?

[bobmyhill]

Would you mind squashing your commits?

Of course! Now squashed, and I added @gassmoeller as co-author on the first commit as he did much of the initial work.

[gassmoeller]

Thanks, I will take a closer look tomorrow. My only concern so far is what Timo pointed out here. These lines are almost what we do in the call to sim.compute_initial_velocity_boundary_constraints(constraints_v) and sim.compute_current_velocity_boundary_constraints(constraints_v), except that the simulator uses the general DoFHandler, while the special lines were using the DoFHandler for the GMG solver that only includes the velocity variable. Does it matter? I dont know immediately, but would want to check if something changes if I put the lines back in.

[bobmyhill]

I ... would want to check if something changes if I put the lines back in.

The code breaks with a "cycle detected in constraints" deal.II error (or it did when I tried it). I also don't know whether this broke something subtle.

[bobmyhill]

@gassmoeller I restored the file in the new commit above to test it.

[gassmoeller]

I agree with the general direction and discussion of this PR, but I have questions on the lines I marked up.

Overall I think:

• you got the fix right
• we need to take care to only disable the manifold on boundaries that are actually deforming (for all solvers!), otherwise we loose a lot of accuracy as @tiannh7 pointed out
• I think there is at least one change in the mesh deformation system that I dont think is necessary see my comment there
• I think there is something unclear happening in the lines I commented on yesterday and that Bob reverted in the latest commit. constraints_v should only contain constraints for the velocity degrees of freedom as numbered by the DoFHandler dof_handler_v, yet we hand over constraints_v to the functions compute_initial_velocity_boundary_constraints and compute_current_velocity_boundary_constraints which use our general DoFHandler to put some constraints in. It doesnt seem to matter as tests pass either way, but it feels wrong to me. @tjhei do you agree, or am I missing something? I would suggest to leave out the changes you reverted Bob (so just stick to the version we have now without removing the weird lines) and we figure it out separately from this PR. It doesnt seem to matter for the tests either way.

[gassmoeller]

I am not sure about this change. These no normal flux constraints are for the multigrid level constraints of mesh deformation of boundaries that allow only tangential mesh movements. Shouldnt we be able to use the manifold for those boundaries as well?

[gassmoeller]

doesnt this also need a special treatment for boundaries with/without mesh deformation? using the manifold for boundaries without mesh deformation and not using the manifold for boundaries with mesh deformation?

[gassmoeller]

same question as for the global coarsening gmg, shouldnt this be split into one call for the tangential velocity boundaries with mesh deformation, and one for the ones without?

[gassmoeller]

and here the same, this is the same operation as in line 1625, just for the level dofs instead of the global dofs. Therefore it should also be split into the boundaries with/without deformation.

[gassmoeller]

this seems like a large change, and one I would like to understand better. But probably only worth looking into when my other comments have been addressed (as this is a gmg test).

[bobmyhill]

Hi @gassmoeller , thanks for the review. Do you want to take this one over? I'm going to be extremely busy for the next two-four weeks, and probably won't have time to pay attention to this before w/c 23rd February.

[gassmoeller]

Hi, sure. I implemented most of my comments and pushed them to this branch. I need to take another look tomorrow if I got everything (GMG level constraints and diffusion are not checked yet).

@alarshi I just remembered I wanted to ask you about this PR. What the PR does is to always align the tangential velocity boundary conditions with the surface, even if it is deformed. I feel like we discussed for your models at some point that that is actually not what we want (we wanted material to be able to move sideways even at a mountain tip in topography, so it would move into/out of the model domain). Is this still something you need? Because then we should talk how to make that compatible with this PR.

[gassmoeller]

I found the remaining difference between AMG and GMG! The tangential boundary conditions for GMG are not set by compute_current_constraints but in StokesMatrixFreeHandlerLocalSmoothingImplementation<dim, velocity_degree>::setup_dofs(), which was not called after mesh deformation. If I call it after mesh deformation (see my last commit) then AMG (top panel) and GMG (bottom panel) deliver identical results:

Things still left to do:

• split the changes in diffusion.cc into a separate PR, because they are just renaming variables
• fix/test this for spherical geometries (I am not yet distinguishing between boundaries with/without manifold in the GMG solver)
• Check the global coarsening GMG

[bobmyhill]

Amazing, thank you @gassmoeller !

[bobmyhill]

By GitHub rules I can't approve what is nominally my PR, but FWIW, I approve of this PR when the testers pass!

Thank you, Rene :)

[gassmoeller]

I updated the changelog entry in the second commit, since it wasnt just the manifold issues that was fixed, it was also that we did not update the constraints during the model run.

Yes this should be ready if someone with the right permissions comes along and pushes the button 😄

[tjhei]

"This branch has conflicts that must be resolved"

[gassmoeller]

Yes that happened in the meantime. I rebased and fixed the conflict but still need to update the tests. Will leave a comment when this is ready.