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 On Thu, 21 Mar 2019 at 05:30, Konrad <[email protected]> wrote: > 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 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 > [email protected]. > For more options, visit https://groups.google.com/d/optout. > -- 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.
