Inhomogeneous neumann bc for step-20

[Jane Lee]

Dear all,

I am trying to solve the equations in step-20 with inhomogeneous neumann bcs on one of the boundaries and getting something very bizarre. 

I have a rectangular domain with the following:

1. Top boundary has homogeneous conditions: this is applied into the weak form like in step-20

2. Side boundaries have no flux conditions on u and I am using DoFTools::make_zero_boundary_constraints for this.

3. Bottom boundary (where the problem is), I need to impose a non-zero value for the neumann condition. 

Some notes at this point:

1. I have done some serious debugging already. 

2. I know the problem is in how I am imposing the bottom boundary condition - I know this because I have done some verification on the rest of the code using method of manufacture solutions (MMS) for the top dirichlet conditions and bottom AND side boundaries having all homogeneous neumann bcs (so not non-zero for the bottom like I need). I know that without the inhomogeneous neumann bc, the code works well with errors as expected when compared to the test analytical solution in the MMS.

3. I am aware that for the condition I am trying to impose on the bottom, u.n = g, where g is a non-zero constant,  the degrees of freedom on the faces are exactly the normal velocities, so I do the following:

I have tried to impose the inhomogeneous neumann bc from the ideas here: https://groups.google.com/forum/#!searchin/dealii/step-20$20neumann|sort:date/dealii/b7WG5noR7FY/z7CLuQBae0MJ

and used something like

std::map<types::global_dof_index, double> boundary_values; 
{
        types::global_dof_index n_dofs = dof_handler.n_dofs();
        std::vector<bool> componentVector(dim + 1, true);
        componentVector[dim] = false;
        std::vector<bool> selected_dofs(n_dofs);
        std::set< types::boundary_id > boundary_ids;
        boundary_ids.insert(2);
        boundary_ids.insert(3);
        DoFTools::extract_boundary_dofs(dof_handler, ComponentMask(componentVector),
                selected_dofs, boundary_ids);

        for (types::global_dof_index i = 0; i < n_dofs; i++) {
            if (selected_dofs[i]) boundary_values[i] = -1.0;
        }
}
MatrixTools::apply_boundary_values(boundary_values,
        system_matrix, solution, system_rhs);

Please refer to the attached file, where you can see the non-constant boundary solution. 

Can anyone explain why this is happening? and what I am doing wrong here?

[Loylick]

ingomogeneous neumann BC are imposed through including them into a stiffness matrix.

You should not use MatrixTools::apply_boundary_values for this purpose I guess.

Look the Step-7 trutorial for the info on how to calculate boundary integral for neumann BC.

[Wolfgang Bangerth]

On 08/28/2018 08:07 AM, Jane Lee wrote:
>
> I am trying to solve the equations in step-20 with inhomogeneous neumann bcs
> on one of the boundaries and getting something very bizarre.

step-20 uses a mixed formulation in which both the pressure and the velocity
(in essence, the gradient of the pressure) are primary variables of the problem.

In cases like this, it is often unclear what exactly you mean when you say
"inhomogenous Neumann boundary conditions". Can you clarify what exactly you
mean there? Do you want to prescribe the velocity, the pressure, or the normal
derivative of the pressure (which equals the normal component of the velocity)?

Best
W.

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

[Loylick]

if you formulate a weak form for your problem you will see a boundary integral from grad(u),

in step-20 it is equal to zero, since you changed the BC you should take into account this

term in a weak form too.

[Jane Lee]

I believe the Neumann conditions are strongly imposed. 

And yes - I realised that inhomogeneous Neumann bc is ambiguous phrasing. 

I mean that I have a conditions k grad p.n =g, or u.n = g equivalently, I think this is in point 3 in my notes in my original post but i think it was unclear what u was. I want to impose a nonzero condition on the normal derivative of the pressure. 

Many thanks!

[Konrad]

Hi Jane Lee,

I recently came across a similar problem.

On Thursday, August 30, 2018 at 12:35:47 PM UTC+2, Jane Lee wrote:

I believe the Neumann conditions are strongly imposed. 

And yes - I realised that inhomogeneous Neumann bc is ambiguous phrasing. 

I mean that I have a conditions k grad p.n =g, or u.n = g equivalently, I think this is in point 3 in my notes in my original post but i think it was unclear what u was. I want to impose a nonzero condition on the normal derivative of the pressure. 

Yes,  Neumann condition is a bit ambiguous here. It is often better to talk about "essential boundary conditions". Essential means that the boundary conditions need to be built in the function space. If you test with v and do an integration by parts on the (non-mixed) problem -div (K*grad p) = f you will see that the normal derivative of p along the boundary appears in the variational formulation. A condition such as K*(grad p) *n=g is then called "natural". A condition on p itself needs to be hard coded into the function space you are solving in and is called "essential".

Now, that you have a mixed form of the equation you do quite the same: you test with a pair of test functions (tau and q) in suitable spaces (step-20: H(div) x L^2) and then you integrate by parts. What you then see is that a condition on p appears in the variational form (hence here it is natural) whereas any condition on n*u~K*(grad p) *n needs to be hard coded into the space H(div) (i.e., it is essential).

Now, if you use like in step-20 Raviart-Thomas elements of lowest order (RT0) you need to keep in mind that the degrees of freedom are edge moments (edge averages in case of RT0). It is possible to set them manually by building them into a constraint matrix like that:

    ConstraintMatrix           constraints;

    constraints.clear ();

    DoFTools::make_hanging_node_constraints (dof_handler, constraints);

    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
    const unsigned int   dofs_per_face   = fe.dofs_per_face;
    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
    std::vector<types::global_dof_index> local_dof_face_indices (dofs_per_face);

   

    typename DoFHandler<dim>::active_cell_iterator
                                    cell = dof_handler.begin_active(),
                                    endc = dof_handler.end();
    for (; cell!=endc; ++cell)
    {
        cell->get_dof_indices(local_dof_indices);

        for (unsigned int face_n=0;
                face_n<GeometryInfo<dim>::faces_per_cell;
                ++face_n)
        {
            cell->face(face_n)->get_dof_indices(local_dof_face_indices);
            if (cell->at_boundary(face_n))
            {
                if (cell->face(face_n)->boundary_id()==id_bottom_boundary)
                {
                    for (unsigned int i = 0; i<dofs_per_face; ++i)
                    {

                        // see documentation for these functions
                        constraints.add_line (local_dof_face_indices.at(i));
                        constraints.set_inhomogeneity (local_dof_face_indices.at(i),
                                value_of_dof);
                    }
                }
                else // add some constraints on other boundaries. Here: zero-flux
                {
                    for (unsigned int i = 0; i<dofs_per_face; ++i)
                    {
                        constraints.add_line (local_dof_face_indices.at(i));
                    }
                }
            }
        }
        ++cell_n;
    }

    constraints.close();

Once you use the DoFTools::make_sparsity_pattern function you need to make sure to keep the inhomogenities:

        DoFTools::make_sparsity_pattern(dof_handler,
                                              dsp,
                                              constraints,
                                              /*keep_constrained_dofs = */ true);

The same in the assembly routine:

            cell->get_dof_indices(local_dof_indices);
            constraints.distribute_local_to_global(local_matrix,
                                                                local_rhs_local,
                                                                local_dof_indices,
                                                                system_matrix,
                                                                system_rhs,
                                                                /* use inhomogeneities for rhs */ true);

This is how I did it. Hope that helps.

[Wolfgang Bangerth]

On 08/30/2018 04:35 AM, Jane Lee wrote:
> I believe the Neumann conditions are strongly imposed.
>
> And yes - I realised that inhomogeneous Neumann bc is ambiguous phrasing.
>
> I mean that I have a conditions k grad p.n =g, or u.n = g equivalently,
> I think this is in point 3 in my notes in my original post but i think
> it was unclear what u was. I want to impose a nonzero condition on the
> normal derivative of the pressure.

Jane: Konrad already gave the correct answer. There is also the function
VectorTools::project_boundary_values_div_conforming:
https://www.dealii.org/9.0.0/doxygen/deal.II/group__constraints.html#ga86427a4e8b3a7e580dabb4e473005288

[Konrad]

Hi Jane Lee and Wolfgang,

If I understand this function correctly it is a function that projects a H(div)-vector field onto the Raviart-Thomas space. While this function is clearly useful it is in my opinion less well-suited to prescribe (scalar) normal fluxes. If you were to do so you would explicitly need knowledge of the normal direction to first construct a vector field which then projected on the Raviart-Thomas space would give you the desired normal trace (same would be true for traction boundary conditions in which case you would need to know a suitable stress tensor).

On Thursday, August 30, 2018 at 6:41:51 PM UTC+2, Wolfgang Bangerth wrote:

Jane: Konrad already gave the correct answer. There is also the function
VectorTools::project_boundary_values_div_conforming:
https://www.dealii.org/9.0.0/doxygen/deal.II/group__constraints.html#ga86427a4e8b3a7e580dabb4e473005288

Best
  W.

Jane, I believe you could do what Wolfgang suggested using VectorTools::project_boundary_values_div_conforming but then you need to construct the vector field first and then project it on the boundary in RTk.

Best,
Konrad

[Wolfgang Bangerth]

> If I understand this function correctly it is a function that projects a
> H(div)-vector field onto the Raviart-Thomas space. While this function is
> clearly useful it is in my opinion less well-suited to prescribe (scalar)
> normal fluxes.

No, that's exactly what it's supposed to do

> If you were to do so you would explicitly need knowledge of the normal
> direction to first construct a vector field which then projected on the
> Raviart-Thomas space would give you the desired normal trace (same would be
> true for traction boundary conditions in which case you would need to know a
> suitable stress tensor).

That is true and one of the limitations of this function. Do you have a
complicated domain where a priori knowing the normal vector is difficult?

[mrjonm...@gmail.com]

I've been trying to use VectorTools::project_boundary_values_div_conforming in a similar model to what Jane Lee is working on. I get a seg fault however when I run it on multiple processes with a parallel::distributed::Triangulation. Is it meant to work for distributed triangulations? 

I don't see any if (cell->is_locally_owned()) checks in the implementation, so I'm suspecting not. 

Thanks,

Jonathan

[Wolfgang Bangerth]

On 08/31/2018 06:06 PM, mrjonm...@gmail.com wrote:
> I've been trying to use VectorTools::project_boundary_values_div_conforming in
> a similar model to what Jane Lee is working on. I get a seg fault however when
> I run it on multiple processes with a parallel::distributed::Triangulation. Is
> it meant to work for distributed triangulations?
>

> I don't see any if(cell->is_locally_owned())checks in the implementation, so
> I'm suspecting not.

Yes, I think that is true. But because the function works one face at a time,
that should not actually be very difficult to add. Would you like to give this
a go and see how far you get? We'd be quite happy to assist with whatever
question you have!

[Jonathan Matthews]

Will do, thanks!

                           www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
--- You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jane Lee]

Hello all, 

Thank you for your messages - I've had some time to run some test cases now. 

With what you are saying with the normal vector, having looked at the documentation, I'm unsure as to why I need to know this? It says that "i.e. the normal components of the solution u and a given f shall coincide. The function f is given by boundary_function " so I have fed in the function f which is a vector field, without explicitly stating the normal (as I assumed it thought the normals were calculated somehow). 

Doing this, I have something like VectorTools::project_boundary_values_div_conforming (dof_handler, 0, FluxFunction<dim>(), 2, constraints), where the 2 is the boundary id. This seems to run fine though I have something funny going on in the corner where the Dirichlet boundary and side boundary (flux condition) meet. Is this because of some order I've messed up when I am imposing the conditions?? (please see photo for the odd instability in the corner of my domain). The solution inside the domain is correct, and this still happens however much I refine the grid. 

For more info, the VectorTools flux conditions are in my setup_dofs code, within costraints.clear()... constraints.close(). After assembling the system, I don't use the distribute local to global function but I instead distribute the constraints onto the solution after I have solved the system. Is there something wrong with the order?? I am getting a similar issue with equations like step-22 but I will post separately as it is different. 

Many thanks

[Wolfgang Bangerth]

> With what you are saying with the normal vector, having looked at the
> documentation, I'm unsure as to why I need to know this? It says that "i.e.
> the normal components of the solution u and a given f shall coincide. The
> function f is given by |boundary_function| " so I have fed in the function f
> which is a vector field, without explicitly stating the normal (as I assumed
> it thought the normals were calculated somehow).

Yes, that works.

What we were trying to say is this: the boundary conditions are of the form
u.n = g
where g is a scalar function -- say, you have a flow of 3 cubic meters per
second per (square meter of surface area). So, conceptually, we should have a
function that takes a scalar function as argument.

But, as a historical accident, the function takes a vector field, say \vec h,
and then internally computes g=\vec h . n. What this implies is that if you
were given the physical quantity g, and the function requires you to provide
an h, then you need to find a way to compute an h out of the g, and the way to
do this is clearly to set
h = g n
but in order to do that, you will need to know n which is a bit cumbersome to
obtain when all your function h is given is the point x at which to evaluate.

On the other hand, oftentimes you are interested in testing that you converge
to the correct solution, and in that case you will simply want to set h to the
exact solution. (This is, indeed, how the "historical accident" happened.)

> Doing this, I have something like
> VectorTools::project_boundary_values_div_conforming (dof_handler, 0,
> FluxFunction<dim>(), 2, constraints), where the 2 is the boundary id. This
> seems to run fine though I have something funny going on in the corner where
> the Dirichlet boundary and side boundary (flux condition) meet. Is this
> because of some order I've messed up when I am imposing the conditions??
> (please see photo for the odd instability in the corner of my domain). The
> solution inside the domain is correct, and this still happens however much I
> refine the grid.
>
>
> For more info, the VectorTools flux conditions are in my setup_dofs code,
> within costraints.clear()... constraints.close(). After assembling the system,
> I don't use the distribute local to global function but I instead distribute
> the constraints onto the solution after I have solved the system. Is there
> something wrong with the order?? I am getting a similar issue with equations
> like step-22 but I will post separately as it is different.

I'm not sure what is happening there. Does the order in which you put the
Dirichlet and the flux conditions into the constraints object matter?

[Jane Lee]

k great - good to know that I am implementing it correctly. 

I'm not sure - I didn't think it would have mattered, but something seems to be going when you also have a Dirichlet condition application (My code works fine with Dirichlet all around the boundary). 

The values I am using are correct, as I am feeding in the exact solution for u.n conditions on the sides and bottom (since I know u) and the exact value also for p. I've of course modified the right hand side to account for extra source terms in the method of manufactured solutions. I believe I have enough boundary conditions - Dirichlet at top and flux condition on the sides and bottom

I do VectorTools::project_boundary_values_div_conforming (dof_handler, 0, ExactFunction<dim>(), 2, constraints) where 2 is the boundary id for sides and bottom. 

For boundary id 1 which is the top, I put in the Dirichlet condition in the weak form:

local_rhs(i) += -(fe_face_values[velocities].value (i, q) *

fe_face_values.normalvector(q) *

boundary_values[q] *

fe_face_values.JxW(q));

I believe that's all you need to do. I wondered whether the value where the side and top boundary meet as shown in that screenshot is getting messed up once I apply the constraints over the solution after I have solved for it?

[Daniel Arndt]

Jane,

 

For more info, the VectorTools flux conditions are in my setup_dofs code, within costraints.clear()... constraints.close(). After assembling the system, I don't use the distribute local to global function but I instead distribute the constraints onto the solution after I have solved the system. Is there something wrong with the order?? I am getting a similar issue with equations like step-22 but I will post separately as it is different.

That sounds a bit strange to me. If you are using a ConstraintMatrix object (or an AffineConstraints object), you should normally take the constraints into account while assembling, i.e. use

constraints.distrbiute_local_to_glocal(cell_matrix, cell_rhs, dof_indices, system_matrix, system_rhs);

and again after solving using

constraints.distribute(solution);

Best,
Daniel

[Jane Lee]

Hi Daniel,

I did indeed have this in. 

I seem to have fixed it when I rewrote everything from scratch... seems sometimes it is quicker to start from the beginning. 

Thanks to all that helped - I have learned a great deal from this post alone and I hope I can contribute as much.