[deal.II] 1D Diffusion Equation

[Pietro Maximoff]

Hello

I modified step-23 to solve this particular diffusion equation:

mu*du/dt  - d/dx(k*du/dx)  + beta*u = f                                        
0 < x < 100

subject to the following BCs:

u(x,0) = 20
q(0) = 0.1 and
T(100) = 0

While the general shape of the resulting plots are correct, some show some 
local minima and 'kinks'. so, I must be doing something wrong. Could you please 
help examine the code [setup_system() and run()] below and point out my errors.

Thanks

Pietro


===
template <int dim>
void DiffusionEquation<dim>::setup_system ()
{
  GridGenerator::hyper_cube (triangulation, 0, 100);
  triangulation.refine_global (4);
  
  std::cout << "Number of active cells: "
        << triangulation.n_active_cells()
            << std::endl;

  dof_handler.distribute_dofs (fe);

  std::cout << "Number of degrees of freedom: "
        << dof_handler.n_dofs()
        << std::endl
        << std::endl;

  sparsity_pattern.reinit (dof_handler.n_dofs(),  dof_handler.n_dofs(), 
dof_handler.max_couplings_between_dofs());
  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
  sparsity_pattern.compress();

  system_matrix.reinit (sparsity_pattern);
  mass_matrix.reinit (sparsity_pattern);
  mass_matrix_zero.reinit (sparsity_pattern);
  laplace_matrix.reinit (sparsity_pattern);
  total_laplace_matrix.reinit (sparsity_pattern);
  
  DiffusionCoefficient<dim> dC;
  MassMatrixCoefficient<dim> mC;
  MassMatrixZeroCoefficient<dim> mzC;
  
  MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3), mass_matrix, 
&mC);
  MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3), 
mass_matrix_zero, &mzC);
  MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3), 
total_laplace_matrix, &dC);
  
  total_laplace_matrix.add(1, mass_matrix_zero);
  
  system_matrix.copy_from (mass_matrix);
  system_matrix.add (theta * time_step, laplace_matrix);

  solution_u.reinit (dof_handler.n_dofs());
  old_solution_u.reinit (dof_handler.n_dofs());
  system_rhs.reinit (dof_handler.n_dofs());

  constraints.close ();
}

template <int dim>
void DiffusionEquation<dim>::run () 
{
  setup_system();

  VectorTools::project (dof_handler, constraints, QGauss<dim>(3), 
InitialValuesU<dim>(), old_solution_u);

  Vector<double> tmp (solution_u.size());
  Vector<double> forcing_terms (solution_u.size());
  
  for (timestep_number=1, time=time_step;
       time<=50;
       time+=time_step, ++timestep_number)
    {
      std::cout << "Time step " << timestep_number
        << " at t=" << time
        << std::endl;
      
      mass_matrix.vmult (system_rhs, old_solution_u);

      total_laplace_matrix.vmult (tmp, old_solution_u);
      system_rhs.add (-(1-theta) * time_step, tmp);

      RightHandSide<dim> rhs_function;
      rhs_function.set_time (time);
      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
                       rhs_function, tmp);
      forcing_terms = tmp;
      forcing_terms *= theta;
      
      rhs_function.set_time (time-time_step);
      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
                       rhs_function, tmp);

      forcing_terms.add ((1-theta), tmp);

      system_rhs.add (time_step, forcing_terms);

      {
        
         QGauss<dim>  quadrature_formula(3);
         FEValues<dim> fe_values (fe, quadrature_formula,
                update_values   | update_gradients |
                update_quadrature_points | update_JxW_values);
            
        Vector<double>       cell_rhs (fe.dofs_per_cell);
        std::vector<unsigned int> local_dof_indices (fe.dofs_per_cell);
        typename DoFHandler<dim>::active_cell_iterator cell = 
dof_handler.begin_active(),
        endc = dof_handler.end();
       
        for (; cell!=endc; ++cell)
        {
            fe_values.reinit (cell);
            cell_rhs = 0;
                //APPLY NEUMANN BOUNDARY CONDITION
            if (cell->at_boundary(0) == true) {
                const double neumann_value = 0.1;
                cell_rhs(0) += neumann_value;
            }
        
            cell->get_dof_indices (local_dof_indices);
            for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
            {
                system_rhs(local_dof_indices[i]) += cell_rhs(i);
            }
        }
    
        BoundaryValuesU<dim> boundary_values_u_function;
        boundary_values_u_function.set_time (time);
      
        std::map<unsigned int,double> boundary_values;
        VectorTools::interpolate_boundary_values (dof_handler,
                          1,
                          boundary_values_u_function,
                          boundary_values);
        MatrixTools::apply_boundary_values (boundary_values,
                        system_matrix,
                        solution_u,
                        system_rhs);
      }
      
      solve_u ();

      total_laplace_matrix.vmult (system_rhs, solution_u);
      system_rhs *= -theta * time_step;

      total_laplace_matrix.vmult (tmp, old_solution_u);
      system_rhs.add (-time_step * (1-theta), tmp);

      system_rhs += forcing_terms;

      
      output_results ();
      
      old_solution_u = solution_u;
    }
}



      
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii