AFC stabilization method

[Shahin Heydari]

Dear all

Hello, I hope you all are safe and healthy.

I am trying to implement algebraic flux correction stabilization on the evolutionary convection-diffusion-reaction in convection-dominated case. In part of the method I need to compute the raw-flux by going through all the nodes of the global matrix one by one. The raw-flux is

     fij = mij * (ul(n,i) - ul(n,j)) - mij * (u(n-1,i) - u(n-1,j))
         - theta * time_step * dij * (ul(n,i) - ul(n,j))
- (1 - theta) * time_step * dij * (u(n-1,i) - u(n-1,j))

ul(n,i): predictor of the solution value u at time step n and in node i

u(n-1,i): the solution u at the previous time step (n-1) in node i

mij: value of mass_matrix at node(i,j)

dij: value of artificial_diffusion at node (i,j)

and then compute the flux-correction alpha_{i,j} and add the f = sum alpha_{i,j} * fij , i,j =1,...,N to the algebraic equation to control the unwanted diffusion. I have computed fij as follow at time step n:

raw_flux_matrix(i,j) =
     mass_matrix(i,j) * (solution_ul(i) - solution_ul(j))
     -
     mass_matrix(i,j) * (old_solution_u(i) - old_solution_u(j))
     -
     time_step * theta * matrix_D(i,j) * (solution_ul(i) - solution_ul(j))
     -
     time_step * (1-theta) * matrix_D(i,j) * (old_solution_u(i) - old_solution_u(j));

but the obtained solution is very smeared and the solution does not do one complete rotation as it supposed to do. I think the problem might be that I have made a mistake in computing mij, dij and ul in the global indices. As I know, when we are inside the cell and want to access the value of global matrix or vector we can compute the value at the quadrature point for example as follow:

 std::vector<double> old_solution_u_values(n_q_points);

fe_values.get_function_values(old_solution_u ,  old_solution_u_values);

Now I am just wondering how I can access the entry of the global matrix and vector in this case.

I would be very thankful and highly appreciated if somebody can help me and give me some advice.

In the attachment I have uploaded the code and also the related article that I have written the code according to. I am eagerly looking forward to hearing any advice and feedback.

Kind Regards

Shahin Heydari

[Shahin Heydari]

  void CDR<dim>::assemble_raw_flux()
  {
    raw_flux_matrix = 0.0;
    for(unsigned int i = 0; i < dof_handler.n_dofs(); ++i)      
      {
          for(auto jt = sparsity_pattern.begin(i); jt != sparsity_pattern.end(i); jt++);
          const auto j = jt->column();
          raw_flux_matrix(i,j) =

[Wolfgang Bangerth]

On 12/10/21 03:23, Shahin Heydari wrote:
>
> raw_flux_matrix(i,j) =
>      mass_matrix(i,j) * (solution_ul(i) - solution_ul(j))
>      -
>      mass_matrix(i,j) * (old_solution_u(i) - old_solution_u(j))
>      -
>      time_step * theta * matrix_D(i,j) * (solution_ul(i) - solution_ul(j))
>      -
>      time_step * (1-theta) * matrix_D(i,j) * (old_solution_u(i) -
> old_solution_u(j));
> but the obtained solution is very smeared and the solution does not do
> one complete rotation as it supposed to do. I think the problem might be
> that I have made a mistake in computing mij, dij and ul in the global
> indices.

I suspect that you won't get anyone here who has the interest,
knowledge, and time to go through your code to check out what is the
problem. So it comes down to how you can debug the problem and my
suggestion always is to pick simple cases: For example, choose an
initial condition that is a step function and for which you can compute
everything on a piece of paper, run one time step, and compare what you
get with what you expected.

In any case, the code snippet you show that loops through the matrix
entries is at least conceptually correct, though I will point out that
this...

for(auto jt = sparsity_pattern.begin(i); jt !=
sparsity_pattern.end(i); jt++);
const auto j = jt->column();

...looks funny to me: The semicolon at the end of the for(...); line
means that the loop has an empty loop body. But then the variable 'jt'
would not actually live outside the loop body, and I would have expected
that you get an error message in the second line. In any case, the
following line that computes raw_flux_matrix is not within the loop body
because there are no curly braces, and so does not seem to be executed
for every 'j'. That doesn't seem right, but I infer that the code does
not compile that way, and so it must not be exactly what you have in
your code.

> As I know, when we are inside the cell and want to access the
> value of global matrix or vector we can compute the value at the
> quadrature point for example as follow:
>  std::vector<double> old_solution_u_values(n_q_points);
> fe_values.get_function_values(old_solution_u ,  old_solution_u_values);

I don't get that comment. I thought you were just working with nodal
values of the solution for the AFC. Why do you need the values at
quadrature points?

> Now I am just wondering how I can access the entry of the global matrix
> and vector in this case.

And I don't understand this question because you have already figured
this out.

Best
W.

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

[Thomas WICK]

Dear Shahin, dear Wolfgang,

I had a closer look into the technical details. The implementation with
regard to deal.II functionalities is fine. Only for
the nonlinear AFC method, the algorithm needs to be worked out.

Best regards,

Thomas

--
---
++++--------------------------------------------++++
Prof. Dr. Thomas Wick
Leibniz Universität Hannover (LUH)
Institut für Angewandte Mathematik (IfAM)
Arbeitsgruppe Wissenschaftliches Rechnen (GWR)

Welfengarten 1
30167 Hannover, Germany

Tel.: +49 511 762 3360
Email: thoma...@ifam.uni-hannover.de
www: https://ifam.uni-hannover.de/wick
www: https://thomaswick.org
++++--------------------------------------------++++
---
--



Am 16.12.21 um 21:04 schrieb Wolfgang Bangerth:

[Shahin Heydari]

Dear Prof. Bangerth and Prof.Wick

Thank you so much for your reply and suggestion as always, I really appreciate your time.

As Prof. Wick has mentioned we have checked the code and the problem was related to the nonlinear part.

Best Regards

Shahin

--
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/7bbd3d1a-829e-486d-a9c2-02134ba29b18%40gmail.com.