ABOUT ASSEMBLING RIGHT HAND SIDE

[Xiang Sun]

Hi, I am implementing the backward Euler method to solve a time-dependent problem like u(n+1)= h*f(n+1)+u(n). n is the nth time step. h is the time step length. f is a function of time. I tried two ways to assemble the system right-hand side(rhs). 

One is to form a mass matrix first, then time the old value of u at the node. This way is similar to Step-26. The code is like

mass_matrix.vmult(system_rhs, old_solution);

The other way is to get the value of u at quadrature points by interpolation, then do the integral in cells to get the local rhs, finally assemble the system rhs by the local rhs. This method is from Step-31 and shown in the following code. 

    for (; cell != endc; ++cell) {

      local_rhs = 0;

      fe_values.reinit(cell);

      fe_values.get_function_values(old_solution, old_sol_values);

      for (unsigned int q = 0; q < n_q_points; ++q ){

        for (unsigned int i = 0; i < dofs_per_cell; ++i) {

          const double phi_i = fe_values.shape_value(i, q);

          local_rhs(i) += old_sol_values[q] * phi_i * fe_values.JxW(q);

        }

        

      }

      cell->get_dof_indices(local_dof_indices);

      for (unsigned int i = 0; i < dofs_per_cell; ++i) {

        system_rhs(local_dof_indices[i]) += local_rhs(i)*0;

      }

    }

I thought these two methods should give the same answer, but they can't. I don't know where I made a mistake. Thank you very much. 

[Xiang Sun]

For correction, the right code for the second method is 

    for (; cell != endc; ++cell) {

      local_rhs = 0;

      fe_values.reinit(cell);

      fe_values.get_function_values(old_solution, old_sol_values);

      for (unsigned int q = 0; q < n_q_points; ++q ){

        for (unsigned int i = 0; i < dofs_per_cell; ++i) {

          const double phi_i = fe_values.shape_value(i, q);

          local_rhs(i) += old_sol_values[q] * phi_i * fe_values.JxW(q);

        }

        

      }

      cell->get_dof_indices(local_dof_indices);

      for (unsigned int i = 0; i < dofs_per_cell; ++i) {

        system_rhs(local_dof_indices[i]) += local_rhs(i);

      }

    }

[Wolfgang Bangerth]

On 4/5/20 12:32 AM, Xiang Sun wrote:
> For correction, the right code for the second method is

> [...]

Are you saying that you solved your own problem? Or are you still unclear what
is happening?

Best
W.


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

[Jean-Paul Pelteret]

Hi Xiang,

You’ve also not mentioned whether you’re comparing vectors before or after you’ve distributed the constraints. Perhaps you get the same result if you call constraints.distribute() on both of your vectors?

Best,
Jean-Paul

> --
> 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/77b269f1-b91b-1dc3-2145-b023242ba43c%40colostate.edu.

[孙翔]

Hi, Wolfgang,

I am still unclear about what is happening. That zero is just a typo. Even if I remove the zero, these two methods can not get the same solution. What I did is just using two methods to calculate the right-hand side generated by time discretization. One way is to form a mass matrix then time the old solution, the other way is to form the local right-hand side in element level then assemble it into system right-hand side array.  

Best,

Xiang

[Wolfgang Bangerth]

On 4/7/20 4:53 PM, 孙翔 wrote:
>
> I am still unclear about what is happening. That zero is just a typo. Even if
> I remove the zero, these two methods can not get the same solution. What I did
> is just using two methods to calculate the right-hand side generated by time
> discretization. One way is to form a mass matrix then time the old solution,
> the other way is to form the local right-hand side in element level then
> assemble it into system right-hand side array.

Right, these two approaches should result in the same solution. So the
question is: how do they differ? Do they result in different right hand side
vectors, or in different solutions? If they differ in the solutions, once you
visualize them, is it obvious that the problem is at the boundary? At hanging
nodes? These are useful tools to figure out where the bug may be.

[孙翔]

Hi, Jean-Paul,

Do you mean the Dirichlet boundary condition? I applied the boundary condition after I form the right-hand side vector. I don't think the boundary condition will lead to different solutions because the right-hand side vector has been different before I add the boundary conditions. 

Best,

Xiang

On Tuesday, 7 April 2020 13:27:17 UTC-7, Jean-Paul Pelteret wrote:

Hi Xiang,

You’ve also not mentioned whether you’re comparing vectors before or after you’ve distributed the constraints. Perhaps you get the same result if you call constraints.distribute() on both of your vectors?

Best,
Jean-Paul

> On 07 Apr 2020, at 22:22, Wolfgang Bangerth <bang...@colostate.edu> wrote:
>
> On 4/5/20 12:32 AM, Xiang Sun wrote:
>> For correction, the right code for the second method is
> > [...]
>
> Are you saying that you solved your own problem? Or are you still unclear what is happening?
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth          email:                 bang...@colostate.edu
>                           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 dea...@googlegroups.com.

[Wolfgang Bangerth]

On 4/7/20 4:58 PM, 孙翔 wrote:
>
> Do you mean the Dirichlet boundary condition? I applied the boundary condition
> after I form the right-hand side vector. I don't think the boundary condition
> will lead to different solutions because the right-hand side vector has been
> different before I add the boundary conditions.

"I don't think" is the wrong approach :-) You are seeing something you are not
expecting. In other words, what you "expect" is not happening, and so you need
to let go of what you think is happening -- you need to *verify* that what you
think is happening is happening! That's the whole mind set of debugging
problems :-)

Best
WB

[孙翔]

Hi, Wolfgang,

Thanks a lot. I didn't apply refinement in my code, so there is no hanging node problem. As you said, there is a big difference around the boundary. I set the boundary value always equal to the initial value which is equal to 1. The M*old_solution method gives the correct answer which is 1, but the second method can not. The attached file is the solution. 

Best,

Xiang

[孙翔]

Hi, Wolfgang,

I applied a Dirichlet boundary condition after I formed system matrixes and vectors, then use the method you taught in the lecture to change the entries in the system matrix and rhs vector. Could I upload the code here in order to make it clear? Thank you very much.

Best,

Xiang 

On Tuesday, 7 April 2020 16:03:10 UTC-7, Wolfgang Bangerth wrote:

[Wolfgang Bangerth]

On 4/7/20 5:31 PM, 孙翔 wrote:
>
> I applied a Dirichlet boundary condition after I formed system matrixes and
> vectors, then use the method you taught in the lecture to change the entries
> in the system matrix and rhs vector. Could I upload the code here in order to
> make it clear?

You can, but you will have to learn to debug these problems yourself.

If the boundary values are wrong, is this the case if the correct boundary
values are nonzero? Does it also happen if the correct boundary are zero? If
it works for zero boundary values, then you should carefully read step-26
again: You can only use apply_boundary_values() once for a given matrix
because it changes the matrix.

[孙翔]

Dear Wolfgang,

I tried the zero boundary and it can get a correct answer. I truly used the function for the matrix once in a time loop. In the new time loop, I have updated the system matrix and vector. I am trying to debug it. Could you help me review the code? If it is a stupid mistake, please just tell me it is and I will try to continue to debug it. If it is due to some function usage or mathematics, please help me with some details because I am just a beginner and really want to learn more from dealii. Thank you very much.

Best,

Xiang

[Wolfgang Bangerth]

On 4/8/20 10:20 AM, 孙翔 wrote:
>
> I tried the zero boundary and it can get a correct answer. I truly used the
> function for the matrix once in a time loop. In the new time loop, I have
> updated the system matrix and vector. I am trying to debug it. Could you help
> me review the code? If it is a stupid mistake, please just tell me it is and I
> will try to continue to debug it. If it is due to some function usage or
> mathematics, please help me with some details because I am just a beginner and
> really want to learn more from dealii.

I can't debug everyone's code here on this mailing list, I'm afraid. I had a
brief look and the setup of the system matrix looks correct/reasonable to me.
But then you describe that you have two methods of which one works and the
other doesn't, and I don't know what this program implements; I also don't
know what the *expected* solution for this program is. So I'm afraid you'll
have to discuss this yourself. But as in the last few days, I'm happy to
provide ideas for how to debug. You'll just have to do it yourself...

Best
W.

[孙翔]

Dear Wolfgang,

Thank you very much. I will try to see what happened by myself.

Best,

Xiang