Hi everyone,
I have made some changes to step26:
1. took theta as 1.
2. multiplied the coefficients cd(volumetric heat capacity ) and kt(heat
transfer coefficient) in front of mass matrix and stiffness matrix to
10^6(J/(m^3·℃)) and 0.34（W/(m·℃)) respectively.
3. imposed Dirichlet boundary conditions (colorize=true, -50 ℃ for boundary
id=0, 0 for boundary id=2,4, and took the original periodic RHS of step-26
to 0)
[image: 微信图片_20240804165104.png]
Then I noticed that positive temperatures appear in the interior of the
domain.
What's the cause of this?
Any help is appreciate.
[image: 微信图片_20240804164229.png]
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1ccb6c79-18d6-406a-a00a-26f5fb993cffn%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, Texas A&M University, 2013
*/
#include <deal.II/base/utilities.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/matrix_tools.h>
#include <fstream>
#include <iostream>
namespace Step26
{
using namespace dealii;
template <int dim>
class HeatEquation
{
public:
HeatEquation();
void run();
private:
void setup_system();
void solve_time_step();
void output_results() const;
void refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level);
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> old_solution;
Vector<double> system_rhs;
double time;
double time_step;
unsigned int timestep_number;
const double theta;
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide()
: Function<dim>()
, period(0.2)
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
private:
const double period;
};
template <int dim>
double RightHandSide<dim>::value(const Point<dim> & p,
const unsigned int component) const
{
(void)component;
AssertIndexRange(component, 1);
Assert(dim == 2, ExcNotImplemented());
const double time = this->get_time();
const double point_within_period =
(time / period - std::floor(time / period));
if ((point_within_period >= 0.0) && (point_within_period <= 0.2))
{
//if ((p[0] > 0.5) && (p[1] > -0.5))
//return 1;
//else
return 0;
}
else if ((point_within_period >= 0.5) && (point_within_period <= 0.7))
{
//if ((p[0] > -0.5) && (p[1] > 0.5))
//return 1;
//else
return 0;
}
else
return 0;
}
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> & /*p*/,
const unsigned int component) const
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return -50;
}
//-=-=-=-=-=-=-=-=-=-=-=- cd and kt
template <int dim>
class Cd : public Function<dim>
{
public:
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
double Cd<dim>::value(const Point<dim> & /*p*/,
const unsigned int component) const
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 1e6;
}
template <int dim>
class Kt : public Function<dim>
{
public:
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
double Kt<dim>::value(const Point<dim> & /*p*/,
const unsigned int component) const
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0.34;
}
//-=-=-=-=-=-=-=-=-=-=-=-
template <int dim>
HeatEquation<dim>::HeatEquation()
: fe(1)
, dof_handler(triangulation)
, time_step(1. / 500)
, theta(1)//(0.5)
{}
template <int dim>
void HeatEquation<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << std::endl
<< "===========================================" << std::endl
<< "Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl
<< std::endl;
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
system_matrix.reinit(sparsity_pattern);
Cd<dim> cd000;
Kt<dim> kt000;
MatrixCreator::create_mass_matrix(dof_handler,
QGauss<dim>(fe.degree + 1),
mass_matrix,
&cd000);
MatrixCreator::create_laplace_matrix(dof_handler,
QGauss<dim>(fe.degree + 1),
laplace_matrix,
&kt000);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void HeatEquation<dim>::solve_time_step()
{
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
SolverCG<Vector<double>> cg(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
preconditioner.initialize(system_matrix, 1.0);
cg.solve(system_matrix, solution, system_rhs, preconditioner);
constraints.distribute(solution);
std::cout << " " << solver_control.last_step() << " CG iterations."
<< std::endl;
}
template <int dim>
void HeatEquation<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "U");
data_out.build_patches();
data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number));
const std::string filename =
"Step26-output/solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
std::ofstream output(filename);
data_out.write_vtk(output);
}
template <int dim>
void HeatEquation<dim>::refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
estimated_error_per_cell,
0.6,
0.4);
if (triangulation.n_levels() > max_grid_level)
for (const auto &cell :
triangulation.active_cell_iterators_on_level(max_grid_level))
cell->clear_refine_flag();
for (const auto &cell :
triangulation.active_cell_iterators_on_level(min_grid_level))
cell->clear_coarsen_flag();
SolutionTransfer<dim> solution_trans(dof_handler);
Vector<double> previous_solution;
previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
triangulation.execute_coarsening_and_refinement();
setup_system();
solution_trans.interpolate(previous_solution, solution);
constraints.distribute(solution);
}
template <int dim>
void HeatEquation<dim>::run()
{
const unsigned int initial_global_refinement = 2;
const unsigned int n_adaptive_pre_refinement_steps = 4;
GridGenerator::hyper_L(triangulation,-1,1,true);
triangulation.refine_global(initial_global_refinement);
setup_system();
unsigned int pre_refinement_step = 0;
Vector<double> tmp;
Vector<double> forcing_terms;
start_time_iteration:
time = 0.0;
timestep_number = 0;
tmp.reinit(solution.size());
forcing_terms.reinit(solution.size());
VectorTools::interpolate(dof_handler,
Functions::ZeroFunction<dim>(),
old_solution);
solution = old_solution;
output_results();
while (time <= 0.5)
{
time += time_step;
++timestep_number;
std::cout << "Time step " << timestep_number << " at t=" << time
<< std::endl;
mass_matrix.vmult(system_rhs, old_solution);
laplace_matrix.vmult(tmp, old_solution);
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>(fe.degree + 1),
rhs_function,
tmp);
forcing_terms = tmp;
forcing_terms *= time_step * theta;
rhs_function.set_time(time - time_step);
VectorTools::create_right_hand_side(dof_handler,
QGauss<dim>(fe.degree + 1),
rhs_function,
tmp);
forcing_terms.add(time_step * (1 - theta), tmp);
system_rhs += forcing_terms;
system_matrix.copy_from(mass_matrix);
system_matrix.add(theta * time_step, laplace_matrix);
constraints.condense(system_matrix, system_rhs);
{
BoundaryValues<dim> boundary_values_function;
Functions::ZeroFunction<dim> f_zero;
boundary_values_function.set_time(time);
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
boundary_values_function,
boundary_values);
VectorTools::interpolate_boundary_values(dof_handler,
2,
f_zero,
boundary_values);
VectorTools::interpolate_boundary_values(dof_handler,
4,
f_zero,
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
solve_time_step();
output_results();
if ((timestep_number == 1) &&
(pre_refinement_step < n_adaptive_pre_refinement_steps))
{
refine_mesh(initial_global_refinement,
initial_global_refinement +
n_adaptive_pre_refinement_steps);
++pre_refinement_step;
tmp.reinit(solution.size());
forcing_terms.reinit(solution.size());
std::cout << std::endl;
goto start_time_iteration;
}
else if ((timestep_number > 0) && (timestep_number % 5 == 0))
{
refine_mesh(initial_global_refinement,
initial_global_refinement +
n_adaptive_pre_refinement_steps);
tmp.reinit(solution.size());
forcing_terms.reinit(solution.size());
}
old_solution = solution;
}
}
} // namespace Step26
int main()
{
try
{
using namespace Step26;
HeatEquation<2> heat_equation_solver;
heat_equation_solver.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
