advection diffusion with periodic boundaries

[Gary Uppal]

Hi,

I'm trying to solve a time dependent advection-diffusion equation with periodic boundary conditions. Just a simple du/dt = D \nabla^2 u - v \dot \grad u for now. I use a large diffusion constant, so stability shouldn't be an issue. The solution behaves normally in the bulk, but some of the mass gets reflected once the mass gets to a boundary. The total mass (integral of u over space) then decreases each time it goes through the boundary in one direction, and increases if the velocity is set to go the other way. Any idea why this might be the case? The periodic boundary works fine when its just diffusion and no advection. 

Thank you,

Gary

[Konrad]

Hi Gary,

I'm trying to solve a time dependent advection-diffusion equation with periodic boundary conditions. Just a simple du/dt = D \nabla^2 u - v \dot \grad u for now. I use a large diffusion constant, so stability shouldn't be an issue. The solution behaves normally in the bulk, but some of the mass gets reflected once the mass gets to a boundary. The total mass (integral of u over space) then decreases each time it goes through the boundary in one direction, and increases if the velocity is set to go the other way. Any idea why this might be the case? The periodic boundary works fine when its just diffusion and no advection. 

Small technical question: What time stepping scheme are you using? You will probably be aware of it but any explicit scheme in this case has time step restrictions for the advective part this will be CFL (note this is not a sufficient condition) which scales like v\delta t / \delta x <= const. For the diffusive part the time step constraint is much more severe since it scales like D\delta t / (\delta x)^2 <= const. I implemented a similar (also periodic scheme) and I did not run into this trouble.

Best,

Konrad 

[Gary Uppal]

Hi Konrad,

I've been using an implicit scheme, but also tried treating the advection term explicitly. Specifically, I'm solving the system

(M_ij + k*d*A_ij + k*B_ij) Uⁿ_j = M_ij U^n-1_j + F^n-1

where M is the mass matrix, A is the laplace matrix, B is the advection matrix, k is the time step, and d is the diffusion constant. F is the initial condition, which I take to be a gaussian at time 0, and then don't include for future time steps. 

The problem seems to be mainly with the boundary condition. Here are images of the initial condition, and the solution after some time. As you see, the field seems to be not passing completely through to the other side; the contours bunch up at the wall edge.

Some more details:

I mainly took the code from the deal.ii tutorial step 26 for the heat equation, and removed the adaptive refinement and added periodic boundary conditions and an advection term. And changed the solver to "bicgstab"

Here's some of my code for adding the periodic boundary condition: 

   constraints.clear ();

  // MAKE CONSTRAINTS

  DoFTools::make_hanging_node_constraints (dof_handler, constraints);

    // ADD PERIODICITY:

      std::vector<GridTools::PeriodicFacePair<typename DoFHandler<2>::cell_iterator> >

      periodicity_vector;

      GridTools::collect_periodic_faces(dof_handler, bid_one, bid_two, direction,

                                        periodicity_vector);

      DoFTools::make_periodicity_constraints<DoFHandler<2> >

      (periodicity_vector, constraints); 

constraints.close ();  

And here's my code for building the advection matrix:

void Simulator::create_advection_matrix(){

  const QGauss<2> quad((fe.degree+1));

  const QGauss<1> face_quad((fe.degree+1));

  FEValues<2> fe_values(fe, quad,

        update_values | update_gradients | update_quadrature_points | update_JxW_values);

  FEFaceValues<2> fe_face_values(fe, face_quad,

        update_values | update_normal_vectors | update_quadrature_points | update_JxW_values);

  const unsigned int dofs_per_cell = fe.dofs_per_cell;

  const unsigned int n_q_points = quad.size();

  const unsigned int n_face_q_points = face_quad.size();

  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);

  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

  typename DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();

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

    cell_matrix = 0;

    fe_values.reinit(cell);

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

      const Tensor<1,2> velVal = advField.value(fe_values.quadrature_point(q_index)); // advection field value 

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

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

          cell_matrix(i,j) += ( fe_values.shape_value(i,q_index) * 

            velVal * fe_values.shape_grad(j, q_index) * 

            fe_values.JxW(q_index) ); 

        }

      } // for cell indicies

      // LOCAL TO GLOBAL:

      cell->get_dof_indices(local_dof_indices);

      constraints.distribute_local_to_global(cell_matrix, local_dof_indices, advection_matrix); 

    } // for number of quadrature points

  } // for each cell

} // Simulator::create_advection_matrix()

Sorry if that's too many details, I hope its clear. Thanks again for your help!

Best,

Gary

--
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/sCzuH711cMQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Konrad]

Hi Gary,

I've been using an implicit scheme, but also tried treating the advection term explicitly. Specifically, I'm solving the system

(M_ij + k*d*A_ij + k*B_ij) Uⁿ_j = M_ij U^n-1_j + F^n-1

where M is the mass matrix, A is the laplace matrix, B is the advection matrix, k is the time step, and d is the diffusion constant. F is the initial condition, which I take to be a gaussian at time 0, and then don't include for future time steps. 

I believe F is a forcing term (that can also contain boundary conditions) and that it should also be multiplied by k, the initial condition should be $U^0$. Is the forcing zero? It seems ok to me although I would suggest to use an IMEX scheme (something like BDF for the advective part and then an implicit formula for the diffusion such as second order Adams-Bashforth combined with Crank-Nicolson). 

Also, what is the ratio between your advective and your diffusive parts (i.e., the Peclet number) and what happens if you decrease/increase it?

 

The problem seems to be mainly with the boundary condition. Here are images of the initial condition, and the solution after some time. As you see, the field seems to be not passing completely through to the other side; the contours bunch up at the wall edge.

Some more details:

I mainly took the code from the deal.ii tutorial step 26 for the heat equation, and removed the adaptive refinement and added periodic boundary conditions and an advection term. And changed the solver to "bicgstab"

Here's some of my code for adding the periodic boundary condition: 

   constraints.clear ();

  // MAKE CONSTRAINTS

  DoFTools::make_hanging_node_constraints (dof_handler, constraints);

    // ADD PERIODICITY:

      std::vector<GridTools::PeriodicFacePair<typename DoFHandler<2>::cell_iterator> >

      periodicity_vector;

      GridTools::collect_periodic_faces(dof_handler, bid_one, bid_two, direction,

                                        periodicity_vector);

      DoFTools::make_periodicity_constraints<DoFHandler<2> >

      (periodicity_vector, constraints); 

constraints.close ();  

That looks ok to me (the template argument in the call of make_periodicity_constraints is not necessary though) .

The assembly for the advection looks ok. Sorry if that was not particularly helpful. I never encountered such a bug but having said that also note that I never worked with periodic BCs in deal.ii (although I worked with periodic BCs).

Cheers,
Konrad

 

[Gary Uppal]

Hi Konrad,

I figured out my error. I was imposing constraints differently on the advection matrix than from the other matrices. I fixed it by not using the constraints.local_to_global() function and just computing the full matrix and using constraints.condense() on the system matrix after adding them together. 

Yes, the forcing term is zero except at the first time step, but you're right, I should just set the initial condition U_0. I had it this way since I will be including a non-zero point source later on; I just needed to get this first part to work. 

Right now I've been using a constant velocity of 10 and a diffusion constant of 15, but it seems to be working well even for smaller diffusion constants. I will probably try and move to an IMEX scheme next though.

Thanks again for your help!

Gary

--

[Konrad]

Hi Konrad,

I figured out my error. I was imposing constraints differently on the advection matrix than from the other matrices. I fixed it by not using the constraints.local_to_global() function and just computing the full matrix and using constraints.condense() on the system matrix after adding them together. 

Yes, that makes sense. I am glad you solved it :-)

Konrad