Hi Gary,
I've been using an implicit scheme, but also tried treating the advection
> term explicitly. Specifically, I'm solving the system
> (Mij + k*d*Aij + k*Bij) Unj = Mij Un-1j + Fn-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.
> [image: image_before.jpg][image: image_after.jpg]
>
> 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) .
>
>
> 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!
>
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
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
