Boundary conditions on components in with FESystem

[Konrad Simon]

Dear all,

I have a code version in which I am using a FESystem (3d) composed of three elements: FE_Nedelec (w), FE_RaviartThomas (u) and DGQ (p). 

The code compiles fine but gives unreasonable results and I do not see what I am doing wrong or what I forgot.

I need to impose (essential) boundary conditions on the first two components: 

• on w: w x n =0 (n is normal)
• on u: u*n=0 (no normal flux)

I assume that I am doing something wrong in the setup of BCs or the system (I guess the solver is likely ruled since two different solvers give exactly the same nonsense). This is the code snippet that I suspect to be buggy. I hope the names are self explanatory (essentially it is a modification of step-32). Note that dim=3 here.

Any help would be appreciated. 

Best,
Konrad

/////////////////////////////////////////////////////////////

  // System and dof setup

  /////////////////////////////////////////////////////////////

  template <int dim>

  void

  BoussinesqModel<dim>::setup_nse_matrices(

    const std::vector<IndexSet> &nse_partitioning,

    const std::vector<IndexSet> &nse_relevant_partitioning)

  {

    nse_matrix.clear();

    LA::BlockSparsityPattern sp(nse_partitioning,

                                nse_partitioning,

                                nse_relevant_partitioning,

                                this->mpi_communicator);

    Table<2, DoFTools::Coupling> coupling(2 * dim + 1, 2 * dim + 1);

    for (unsigned int c = 0; c < 2 * dim + 1; ++c)

      {

        for (unsigned int d = 0; d < 2 * dim + 1; ++d)

          {

            if (c < dim)

              {

                if (d < 2 * dim - 1)

                  coupling[c][d] = DoFTools::always;

                else

                  coupling[c][d] = DoFTools::none;

              }

            else if ((c >= dim) && (c < 2 * dim))

              {

                coupling[c][d] = DoFTools::always;

              }

            else if (c == 2 * dim)

              {

                if ((d >= dim) && (d < 2 * dim))

                  coupling[c][d] = DoFTools::always;

                else

                  coupling[c][d] = DoFTools::none;

              }

          }

      }

    DoFTools::make_sparsity_pattern(nse_dof_handler,

                                    coupling,

                                    sp,

                                    nse_constraints,

                                    false,

                                    Utilities::MPI::this_mpi_process(

                                      this->mpi_communicator));

    sp.compress();

    nse_matrix.reinit(sp);

  }

  template <int dim>

  void

  BoussinesqModel<dim>::setup_dofs()

  {

    TimerOutput::Scope timing_section(

      this->computing_timer, "BoussinesqModel - setup dofs of systems");

    /*

     * Setup dof handlers for nse and temperature

     */

    std::vector<unsigned int> nse_sub_blocks(3, 0);

    nse_sub_blocks[1] = 1;

    nse_sub_blocks[2] = 2;

    nse_dof_handler.distribute_dofs(nse_fe);

    /*

     * Reduce bandwidth for ILU preconditioning

     */

    DoFRenumbering::Cuthill_McKee(nse_dof_handler);

    DoFRenumbering::block_wise(nse_dof_handler);

    /*

     * Count dofs

     */

    std::vector<types::global_dof_index> nse_dofs_per_block(3);

    DoFTools::count_dofs_per_block(nse_dof_handler,

                                   nse_dofs_per_block,

                                   nse_sub_blocks);

    const unsigned int n_w = nse_dofs_per_block[0], n_u = nse_dofs_per_block[1],

                       n_p = nse_dofs_per_block[2],

                       n_T = temperature_dof_handler.n_dofs();

    /*

     * Comma separated large numbers

     */

    std::locale s = this->pcout.get_stream().getloc();

    this->pcout.get_stream().imbue(std::locale(""));

    /*

     * Print some mesh and dof info

     */

    this->pcout << "Number of active cells: "

                << this->triangulation.n_global_active_cells() << " (on "

                << this->triangulation.n_levels() << " levels)" << std::endl

                << "Number of degrees of freedom: " << n_w + n_u + n_p + n_T

                << " (" << n_w << " + " << n_u << " + " << n_p << " + " << n_T

                << ")" << std::endl

                << std::endl;

    this->pcout.get_stream().imbue(s);

    /*

     * Setup partitioners to store what dofs and matrix entries are stored on

     * the local processor

     */

    {

      nse_index_set = nse_dof_handler.locally_owned_dofs();

      nse_partitioning.push_back(nse_index_set.get_view(0, n_w));

      nse_partitioning.push_back(nse_index_set.get_view(n_w, n_w + n_u));

      nse_partitioning.push_back(

        nse_index_set.get_view(n_w + n_u, n_w + n_u + n_p));

      DoFTools::extract_locally_relevant_dofs(nse_dof_handler,

                                              nse_relevant_set);

      nse_relevant_partitioning.push_back(nse_relevant_set.get_view(0, n_w));

      nse_relevant_partitioning.push_back(

        nse_relevant_set.get_view(n_w, n_w + n_u));

      nse_relevant_partitioning.push_back(

        nse_relevant_set.get_view(n_w + n_u, n_w + n_u + n_p));

    }

    /*

     * Setup constraints and boundary values for NSE. Make sure this is

     * consistent with the initial data.

     */

    {

      nse_constraints.clear();

      nse_constraints.reinit(nse_relevant_set);

      DoFTools::make_hanging_node_constraints(nse_dof_handler, nse_constraints);

      /*

       * Lower boundary (id=0)

       */

      VectorTools::project_boundary_values_curl_conforming_l2(

        nse_dof_handler,

        /*first vector component */ 0,

        Functions::ZeroFunction<dim>(dim),

        /*boundary id*/ 0,

        nse_constraints);

      VectorTools::project_boundary_values_div_conforming(

        nse_dof_handler,

        /*first vector component */

        3,

        Functions::ZeroFunction<dim>(dim),

        /*boundary id*/ 0,

        nse_constraints);

      /*

       * Upper boundary (id=1)

       */

      VectorTools::project_boundary_values_curl_conforming_l2(

        nse_dof_handler,

        /*first vector component */ 0,

        ConstantFunction<dim>(1, dim),

        /*boundary id*/ 1,

        nse_constraints);

      VectorTools::project_boundary_values_div_conforming(

        nse_dof_handler,

        /*first vector component */

        3,

        Functions::ZeroFunction<dim>(dim),

        /*boundary id*/ 1,

        nse_constraints);

      nse_constraints.close();

    }

   

    /*

     * Setup the matrix and vector objects.

     */

    setup_nse_matrices(nse_partitioning, nse_relevant_partitioning);

    nse_rhs.reinit(nse_partitioning,

                   nse_relevant_partitioning,

                   this->mpi_communicator,

                   true);

    nse_solution.reinit(nse_relevant_partitioning, this->mpi_communicator);

    old_nse_solution.reinit(nse_solution);

  }

[Wolfgang Bangerth]

On 11/6/20 10:11 AM, Konrad Simon wrote:
>
> I have a code version in which I am using a FESystem (3d) composed of three
> elements: FE_Nedelec (w), FE_RaviartThomas (u) and DGQ (p).
>
> The code compiles fine but gives unreasonable results and I do not see what I
> am doing wrong or what I forgot.

Instead of trying to divine what the code is doing (and apparently doing
wrong), let me outline the way I would think about this:

* You say that the results are "unreasonable". That can mean a lot of
different things. Try to explain for yourself what you see, and understand
what part looks wrong. Is it the boundary values? That is, does the solution
you see visualized violate the boundary conditions you wanted to impose? If
so, that's clearly an area of concern.

* If the solution does seem to satisfy the boundary conditions, but is wrong
in the interior of the domain, then the problem lies elsewhere.

* In any case, if you can't find where things go wrong, try solving a
decoupled problem where you assemble a matrix in which each variable appears
in only one equation. This helps you identify *which boundary condition*, or
*which equation* is wrongly implemented.

In the end, debugging comes down to generating hypotheses about what could be
wrong, and then testing each hypothesis in detail until you can confirm or
deny that your hypothesis is the cause of the problem. Creating hypotheses
often comes down to carefully looking at the solution and verbalizing for
yourself what you see and how it (qualitatively) differs from what you expect.
In practice, this often leads to faster ways to debugging than randomly poking
at individual pieces of your code.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Konrad Simon]

Dear Wolfgang,

Many thanks for your reply. I am usually doing as you suggest and I was testing many debugging hypothesis during the last two weeks. I actually found two minor bugs but the main problem persists. I do not know any further so I will be a bit more detailed on where and when the issue arises and I will post an example code based on step-20 to reproduce the issue. It might be an issue within deal.ii which I usually do not assume to be the case.

A few words about the tests I performed:

The mathematical problem I am aiming at is a (dynamic) stokes problem on a hyper_shell in 3D with divergence zero constraint. However, in contrast to step-32 I do not use a H1-L2 formulation for velocity-pressure. I use an H(curl)-H(div)-L2 form for vorticity-velocity-pressure. For mathematical reasons I can not use no-slip BCs but I can use (essential) BCs to constrain the tangential curl (zero in my case) and the normal flux (also zero) on all boundaries. I am looking at the test case of two rising warm bubbles (hyper shell) and one bubble (hyper cube).

My boundary conditions rule out rigid rotations as does the modification of step-32 (no-slip on lower boundary of shell). The pressure in my case is determined up to a constant which I catch in the iterative solver with projections onto the pressure space of L^2 with zero mean. I am using an FESystem with Nedelec-RaviartThomas-DGQ elements that are stable (all lowest order). 

A few observations:

1. The problem in the above test appears only on hyper_shell and hyper_ball-like domains in 3D with Nedelec-RaviartThomas-DGQ. On the cube everything looks OK. 

2. The modification of step-32 I use as a sanity check (the BCs are different but the dynamics should at least be similar). There the problem does not arise. Note that there the problem uses stable Taylor-Hood elements while I use Nedelec-RaviartThomas-DGQ.

3. I tried to reproduce the issue with a simpler problem. For this I modified step-20 where we have an exact solution. I am solving the problem using the exact solution in two forms to compare (natural BCs on p and none on u as well as essential BCs on u which means no-flux and none on p). Both solutions should be the same apart form minor inaccuracies. I do so on a cube and on the hyper_shell.

3a. In 2D the problem does not arise at all. Not on the cube neither on the hyper_shell

3b. In 3D the problem arises in both BC cases (natural on p and essential on u) on the hyper_shell. On the cube everything is fine.

3c. I also project the exact solution with Vectortools::project and add it to the vtu-output to compare. Funnily, the projected exact solution on the hyper_shell in 3D (not in 2D) also looks bad.

3d. I attached some plots of that test case.

3e. I attached a code (deal.ii v9.2) to reproduce the issue (modification of step-20). There are some self-explanatory switches in the main function

Question/Hypothesis:

I am using Vectortools::project_boundary_values_div_conforming with the Raviart-Thomas part of the FESystem. This function must rely on consistent face orientations somehow (note that this is not critical in H^1). I assume that this is also true for the L^2-projections on the RaviartThomas space in the Vectortools::project function since this probably assembles a mass matrix and solves a positive definite system. Therefore, mappings are involved and they need consistent orientations. 

Is my intuition on this correct?

Sorry for bugging you again with this. I know it is unlikely that the error is on deal's side but I ran out of ideas. As said, please see the attached code and plots in 2D and 3D.

Best,

Konrad