Make sure to initialize multigrid preconditioner for steady linear structure

[kronbichler]

This addresses another issue I found in the context of the errors reported in #666 and #664, which is that the multigrid preconditioner did not set up its constituents (in particular the diagonals of the Jacobi/Chebyshev smoother) for the case that we used a steady linear elasticity problem as typical for our FSI problems, see e.g.

exadg/applications/fluid_structure_interaction/cylinder_with_flag/application.h

Lines 671 to 679 in 27a6077

	Parameters & param = this->ale_elasticity_param;
	param.problem_type = ProblemType::Steady;
	param.body_force = false;
	param.pull_back_body_force = false;
	param.large_deformation = false;
	param.pull_back_traction = false;
	param.degree = this->param.mapping_degree;

or similarly

exadg/applications/fluid_structure_interaction/perpendicular_flap/application.h

Lines 531 to 537 in 27a6077

	Parameters & param = this->ale_elasticity_param;
	param.problem_type = ProblemType::Steady;
	param.body_force = false;
	param.pull_back_body_force = false;
	param.large_deformation = false;
	param.pull_back_traction = false;

Reflecting to how we set up the Poisson solver,

exadg/include/exadg/poisson/preconditioners/multigrid_preconditioner.cpp

Lines 55 to 62 in 27a6077

	Base::initialize(mg_data,
	grid,
	multigrid_mappings,
	fe,
	data.operator_is_singular,
	dirichlet_bc,
	dirichlet_bc_component_mask,
	true /* initialize_preconditioners */);

we also need to make sure to build everything for the preconditioner choice that does not call the respective update() function. @davidscn @jh66637 it would be great if you could check this and see if this or adaptations of this help in your case.

[davidscn]

I now applied all patches manually, including #667 (and the deal.II fix) as well as the fix here:

Building and running the same case in release works flawless. However, building and running the testcase in debug mode runs into another exception. Similar to #667 I removed now all manifold_ids. Not entirely sure that this is now the way to go. With this, I get:

  Viscous step:
  Solver viscous step:                       CG
  Maximum number of iterations:              1000
  Absolute solver tolerance:                 1.0000e-08
  Relative solver tolerance:                 1.0000e-06
  Maximum size of Krylov space:              30
  Preconditioner viscous step:               PointJacobi
  Update preconditioner viscous:             false

--------------------------------------------------------
An error occurred in line <378> of file <~/dealii/include/deal.II/grid/tria_objects.h> in function
    unsigned int dealii::internal::TriangulationImplementation::TriaObjects::n_objects() const
The violated condition was: 
    ::dealii::deal_II_exceptions::internals::compare_for_equality(cells.size(), manifold_id.size() * 2 * this->structdim)
Additional information: 
    Two sizes or dimensions were supposed to be equal, but aren't. They
    are 1073741824 and 0.

Stacktrace:
-----------
#0  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: dealii::internal::TriangulationImplementation::TriaObjects::n_objects() const
#1  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: dealii::TriaAccessorBase<2, 2, 2>::operator++()
#2  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: dealii::TriaRawIterator<dealii::CellAccessor<2, 2> >::operator++()
#3  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::begin_quad(unsigned int) const
#4  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::begin_active_quad(unsigned int) const
#5  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::begin_active(unsigned int) const
#6  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::active_cell_iterators() const
#7  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::update_reference_cells()
#8  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::reset_policy()
#9  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: dealii::Triangulation<2, 2>::create_triangulation(std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > > const&, std::vector<dealii::CellData<2>, std::allocator<dealii::CellData<2> > > const&, dealii::SubCellData const&)
#10  ~/dealii/build/lib/libdeal_II.g.so.9.6.0-pre: void dealii::GridGenerator::subdivided_hyper_rectangle<2, 2>(dealii::Triangulation<2, 2>&, std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::Point<2, double> const&, dealii::Point<2, double> const&, bool)
#11  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: ExaDG::FluidFSI::Application<2, double>::create_triangulation(dealii::Triangulation<2, 2>&)
#12  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}::operator()(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&) const
#13  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: void std::__invoke_impl<void, ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}&, dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&>(std::__invoke_other, ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}&, dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int&&, std::vector<unsigned int, std::allocator<unsigned int> > const&)
#14  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: std::enable_if<is_invocable_r_v<void, ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}&, dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&>, void>::type std::__invoke_r<void, ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}&, dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&>(ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}&, dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int&&, std::vector<unsigned int, std::allocator<unsigned int> > const&)
#15  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: std::_Function_handler<void (dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&), ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)::{lambda(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)#1}>::_M_invoke(std::_Any_data const&, dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int&&, std::vector<unsigned int, std::allocator<unsigned int> > const&)
#16  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: std::function<void (dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)>::operator()(dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&) const
#17  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: void ExaDG::GridUtilities::create_triangulation<2>(std::shared_ptr<dealii::Triangulation<2, 2> >&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::Triangulation<2, 2>::cell_iterator>, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::Triangulation<2, 2>::cell_iterator> > >&, ompi_communicator_t* const&, ExaDG::GridData const&, bool, std::function<void (dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::Triangulation<2, 2>::cell_iterator>, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::Triangulation<2, 2>::cell_iterator> > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)> const&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)
#18  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: void ExaDG::GridUtilities::create_triangulation_with_multigrid<2>(ExaDG::Grid<2>&, ompi_communicator_t* const&, ExaDG::GridData const&, bool, std::function<void (dealii::Triangulation<2, 2>&, std::vector<dealii::GridTools::PeriodicFacePair<dealii::Triangulation<2, 2>::cell_iterator>, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::Triangulation<2, 2>::cell_iterator> > >&, unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&)> const&, std::vector<unsigned int, std::allocator<unsigned int> >)
#19  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: ExaDG::FluidFSI::Application<2, double>::create_grid(ExaDG::Grid<2>&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)
#20  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: ExaDG::FluidFSI::ApplicationBase<2, double>::setup(std::shared_ptr<ExaDG::Grid<2> >&, std::shared_ptr<dealii::Mapping<2, 2> >&, std::shared_ptr<ExaDG::MultigridMappings<2, double> >&)
#21  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: ExaDG::FSI::SolverFluid<2, double>::setup(std::shared_ptr<ExaDG::FluidFSI::ApplicationBase<2, double> >, ompi_communicator_t*, bool)
#22  ~/exadg/build_debug/applications/fluid_structure_interaction/perpendicular_flap/solver_precice: ExaDG::FSI::preCICE::DriverFluid<2, double>::setup_fluid_and_ale()
--------------------------------------------------------

[davidscn]

Maybe as an additional information: If I don't remove the set_all_manifold_ids function, but replace it

-        tria.set_all_manifold_ids(0);
+        tria.set_all_manifold_ids(dealii::numbers::flat_manifold_id);

the code runs in debug as well as release mode. I'm not sure what is supposed to work and what not. From the user perspective, I would say in case I don't have a sphere or something I want the code to account for, it should work without setting the flat manifold explicitly.

[kronbichler]

From the user perspective, I would say in case I don't have a sphere or something I want the code to account for, it should work without setting the flat manifold explicitly.

I agree, that's what we should be trying to do.

[nfehn]

Thanks for figuring this out! I am not sure if the suggested change is how we should approach this technically (For unsteady linear problems, the present change could lead to the use of "wrong preconditioners" that do not take into account the time derivative term, or with a wrong scaling. I think we can only do a proper setup of the preconditioner prior to solving the problem for the first time). Generally, I thought about this problem in more detail some time ago during PR #654. For now, I think we should address the problem as in #654. @kronbichler do you know whether the problem you address here would be fixed by the code changes of PR #654?

[kronbichler]

I tested it on one of the problems of #664 and it does also work (I had tried that option as well, but forgot about the connection with #654); I agree with using #654 for now.

You are right that the unsteady case would set up the preconditioner in a wrong setting, so the correct setting would be to pass in a check like we had in #654, nonlinear or data.unsteady, which then would repeat that old code (removed in #654).

I did test the code of #654 and it seemed that it did unnecessary updates, but I would need to investigate that further and is maybe a separate issue. Let me do that later and close this for the time being.

[richardschu]

Regarding "and it seemed that it did unnecessary updates, but I would need to investigate that further and is maybe a separate issue": You mean that i) the change proposed in #654 is not the most efficient way to update the data structures, or ii) that it covers also cases, where the update is indeed not needed? In case you meant ii): afaik, we have Unsteady (linear + nonlinear), QuasiStatic (steady, nonlinear), Steady (linear). So the if-condition only affects the Steady case, which is always linear. So it must be i) you mean, right? @kronbichler

[richardschu]

The if I removed was
further simplified, see the comment

[kronbichler]

You mean that i) the change proposed in #654 is not the most efficient way to update the data structures, or ii) that it covers also cases, where the update is indeed not needed? In case you meant ii): afaik, we have Unsteady (linear + nonlinear), QuasiStatic (steady, nonlinear), Steady (linear

As I said, I did not check closely, I was referring to case ii) in the sense when one needs to refresh the preconditioner, but it seems this was only in my particular example but is captured correctly in the FSI example cylinder_with_flag. I need to investigate the conditions upon which the update is forced, but that is, as I said, a different topic and unrelated to #654.

The if I removed was further simplified, [...]

Yes, which is why I wrote "removed in #654".

[nfehn]

Prinicipally, the idea is that the user has full control regarding a repeated update of preconditioners during a nonlinear solution process or during time stepping. So if one knows there is no need to update a preconditioner (because of the PDE model problem), the user should deactivate this in the parameter settings to avoid "unnecessary" preconditioner updates.

The code initializes/updates preconditioners automatically (i.e. independently of user parameters controlling the update of preconditioners as described above) only in the setup-phase. During the setup phase, we should make sure that the code does not perform "unnecessary" updates. I would consider such a behavior with "unnecessary updates" sub-optimal or erroneous.