Help with step-7 Neumann boundary conditions

[Ali Seddiq]

Hi everyone,

I am trying to mimic step-7 implementation of Neumann boundary condition in the case of a vector-valued problem.

I did some modifications as shown in the below, but it doesn't work as expected;

const FEValuesExtractors::Scalar pressure(0);

.....
        cell_rhs(i) +=
          (fe_face_values[pressure].value(i, q_point) * 
           neumann_value *                                             
           fe_face_values.JxW(q_point));          

I tried to replicate step-20 and related documentation on vector problems but haven't yet been able to figure out the solution.

So any help is much appreciated.

Thank you,

Ali

[Wolfgang Bangerth]

Ali -- a good place to always start is by asking what "it doesn't work as
expected" precisely is. We don't know what the problem is you are trying to
solve, we don't know whether you get compiler errors, whether you get runtime
errors, whether the solution looks wrong, in which way it looks wrong, etc.

Debugging typically starts with asking *yourself* precisely what doesn't work,
and what that might mean for where the bug is.

Best
W.

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

[Ali Seddiq]

Dear Wolfgang,

Thanks for a quick reply, and sorry for ambiguity. To be more precise (trying to be) I should say the lines I modified with the aim of creating a so-called "view" do not do that for my intended variable. My question mainly was that if using the [ ] operator (that I added) does the job of identifying the desired field on the faces? Is the usage correct?

I am also a bit confused as fe.system_to_component_index(i) seems to do the same but not for the faces I guess. So not sure which approach should be taken.

Thank you,

Ali

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/01451b66-7c40-3eb3-53c3-3024eaaad3db%40colostate.edu.

[Wolfgang Bangerth]

On 1/31/22 11:41, Ali Seddiq wrote:
>
> Thanks for a quick reply, and sorry for ambiguity. To be more precise
> (trying to be) I should say the lines I modified with the aim of
> creating a so-called "view" do not do that for my intended variable. My
> question mainly was that if using the [ ] operator (that I added) does
> the job of identifying the desired field on the faces? Is the usage correct?

Yes. The call in your example code considers the first vector component
of the i'th shape function at quadrature point q_point.

> I am also a bit confused as fe.system_to_component_index

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassFiniteElement.html%23a86644fe67824373cd51e9ff7fca94f8c&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C66a9d6d095e7444d971c08d9e4e94b94%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637792512914352482%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=%2BiFE%2FIimkawMXkp54Z7Kci%2BA7AYn%2B6EJg%2BFeTE5LVYQ%3D&reserved=0>(i)

> seems to do the same but not for the faces I guess. So not sure which
> approach should be taken.

system_to_component index is a property of the finite element, not of
the FEValues or FEFaceValues class. As such, it can be used as well.
step-8 shows how it is used for FEValues, but the same could be done for
FEFaceValues.

[Ali Seddiq]

Dear Wolfgang,

Sorry to raise this topic again. But when I add the boundary term to the to cell_rhs on assembly_system, I receive the following message, based on which most likely this addition leads to an invertible matrix:

An error occurred in line <457> of file </usr/local/include/deal.II/lac/solver_cg.h> in function
    void dealii::SolverCG<VectorType>::solve(const MatrixType&, VectorType&, const VectorType&, const PreconditionerType&) [with MatrixType = dealii::SparseMatrix<double>; PreconditionerType = dealii::PreconditionIdentity; VectorType = dealii::Vector<double>]
The violated condition was:
    false
Additional information:
Iterative method reported convergence failure in step 25. The residual in the last step was 1.55891e-08.

This error message can indicate that you have simply not allowed a sufficiently large number of iterations for your iterative solver to converge. This often happens when you increase the size of your problem. In such cases, the last residual will likely still be very small, and you can make the error go away by increasing the allowed number of iterations when setting up the SolverControl object that determines the maximal number of iterations you allow.

The other situation where this error may occur is when your matrix is not invertible (e.g., your matrix has a null-space), or if you try to apply the wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric or not positive definite). In these cases, the residual in the last iteration is likely going to be large.

This is while if I change the pressure, naming it as U2 to another scalar field variable , naming it as U1, it works, though I get zero iteration for U1, and in fact it does not solve for U1.

For example if I write:

fe_face_values[U1].value(i, q_point) *...    (U1: second scalar variable)

It leads to the error, but if I write:

fe_face_values[U1].value(i, q_point) *... (U1: first scalar variable)

does  a number of iterations for U1 , and 0 iteration (probably also meaning no solution) for U2.

As a bit of further information, my weak form simply ends up to:

dt(M+A)*U1 + A'* U2 = dt*(F_{U1^n-1}) + F_U2  

U1 and U2 as two scalar unknowns.

And, the boundary term in matrix-vector form supposed to be my F_U2.

I considered the following block matrices and vectors in forming the linear system (and solver):

(M+A): system_matrix.block(0,0)

A': system_matrix.block(1,1)

U1: solution.block(0)

U2: solution.block(1)

F_U1: system_rhs.block(0)

F_U2: system_rhs.block(1)

F_U1 would be integration over the domain with  phi_i_U1* U1^(n_1)

F_U2 would be integration over the boundary with phi_i_U2*neumann_value

cell_rhs(i) += time_step* phi_i_U1 * old_solution_values * fe_values.JxW(q_point)  ///as  (F_U1)

and there after boundary term is added as:

cell_rhs(i) += neumann_value *fe_face_values[U2].value(i, q_point)  * fe_values.JxW(q_point)  /// as (F_U2)

and the solve function which could be one of the problems is as the following:

        {
                 tmp += system_rhs.block(0);

                 SolverControl            solver_control(solution.block (0).size(), 1e-8 * system_rhs.l2_norm());
                 SolverCG<Vector<double>> cg (solver_control);
                 cg.solve(system_matrix.block(0, 0), solution.block (0), tmp, PreconditionIdentity());

                 std::cout << "     " << solver_control.last_step() << " CG iterations for U1." << std::endl;
         }

         {
                 tmp2 += system_rhs.block(1);

                 SolverControl            solver_control2(solution.block (1).size(), 1e-8 * system_rhs.l2_norm());
                 SolverCG<Vector<double>> cg (solver_control2);
                 cg.solve(system_matrix.block(1, 1), solution.block (1), tmp2, PreconditionIdentity());

                 std::cout << "     " << solver_control2.last_step() << " CG iterations for U2." << std::endl;

         }

I did quite a number of changes in the code to see the outcome but it seems I got stuck on an apparently obvious problem and unfortunately haven't been able to figure it out yet.

I am also not sure if simply addition of boundary term to cell_rhs leads to the correct recognition as system_rhs.block(1) in the solve?

Would you or anyone else please help me with this problem?

Thanks very much,

Ali

[Wolfgang Bangerth]

On 2/8/22 09:12, Ali Seddiq wrote:
> This is while if I change the pressure, naming it as U2 to another
> scalar field variable , naming it as U1, it works, though I get zero
> iteration for U1, and in fact it does not solve for U1.
> For example if I write:
>
> fe_face_values[U1].value(i, q_point) *...    (U1: second scalar variable)
>
> It leads to the error, but if I write:
> fe_face_values[U1].value(i, q_point) *... (U1: first scalar variable)
>
> does  a number of iterations for U1 , and 0 iteration (probably also
> meaning no solution) for U2.

Ali -- we don't know what your code looks like (and even if we did, we
would not have the time to debug it), but the paragraph above suggests
to me that you do not actually understand what the various parts of your
code are doing.

My suggestion to you is to make sure you have the simplest possible case
that fails. That means that you find out what the difference between U1
and U2 is, and why one fails and the other does not. In general, I find
that copying the current version of a program to a safe place, and then
slashing everything from it that is not necessary to show the bug is a
valid strategy for debugging. Make the program as small as you can,
since this reduces the number of places where the bug can be.

[Ali Seddiq]

Dear Wolfgang,

Thanks very much for your advice. I definitely should and will follow that. But may I still narrow my question for a clarification (with upcoming complexity), and as a quick question ask if adding the boundary term through cell_rhs +=... (as above, and very similar to step-7 in terms of implementation, with minor modifications for a vector problem) leads it to get recognized as the second (row) component of the BlockVector F (system_rhs.block(1)) at the right hand side of the linear system? I mean are there any other general requirements than adding "fe_face_values[U2]..." to the cell_rhs for this purpose?

Thank you,

Ali

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/a9_hHimuNLQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/af8b61af-423b-7355-e9ae-4269367a074a%40colostate.edu.

[Wolfgang Bangerth]

On 2/11/22 09:54, Ali Seddiq wrote:
>
> Thanks very much for your advice. I definitely should and will follow that.
> But may I still narrow my question for a clarification (with upcoming
> complexity), and as a quick question ask if adding the boundary term through
> cell_rhs +=... (as above, and very similar to step-7 in terms of
> implementation, with minor modifications for a vector problem) leads it to get
> recognized as the second (row) component of the BlockVector F
> (system_rhs.block(1)) at the right hand side of the linear system? I mean are

> there any other general requirements thanadding "fe_face_values[U2]..." to the
> cell_rhs for this purpose?

Ali -- there is really not enough information on your question to answer this
way or that way.

You have to learn how to debug these sorts of things. If I understand what you
are asking correctly, you want to know whether a particular syntax
is/can/should be used to put a specific entry into the right hand side vector.
This is not hard to check. Modify your program so that it has, for example,
only one or at most four cells (in 2d). Modify your program so that the
assembly function only has code for the one term for which you want to find
out what it does. Modify your program so that you output the right hand side
vector right after assembly. Then compare what you get with what you expect to
see. If it is the same, you just proved yourself that your code is doing what
you were hoping. If it is not the same, then it is time to reconsider.

These sorts of debugging skills are something you have to acquire. It takes
time to ask others questions, and it often takes substantial time to obtain
answers. It is usually faster to just find these answers yourself. You also
learn more this way.

[Ali Seddiq]

Dear Wolfgang,

Thank you many times for your advice. In fact, simplifying the problem led me to find the error in a totally unexpected (to me!) place in the code and made me to define my problem more strictly.

Best Regards,

Ali

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/24e70def-df38-ad59-3a0b-efc31abb76fd%40colostate.edu.