In the context of a parallel implementation of step-26, where I inspired
myself in step-32, I am trying to set a different constant in the initial
and boundary conditions that I apply to my system. Below is an image of the
result when I set the initial system to 1 instead of 0. I tried to
investigate the implementation of a rhs. Somehow the conditions appear to
be inhomogeneous and therefore the rhs is modified. The only improvement
that I got was when I muted this lines in my code:
if (constraints.is_inhomogeneously_constrained(
data.local_dof_indices[i]))
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
data.matrix_for_bc(j, i) +=0 ;
/*
(scratch.phi_T[i] * scratch.phi_T[j] *
(use_bdf2_scheme ? ((2 * time_step + old_time_step) /
(time_step + old_time_step)) :
1.) ) *
scratch.fe_values.JxW(q);
*/
}
but I am still not able to set the entire system initial temperature to 1
or any other initial value.
Can you help me?
Thank you!
[image: Screenshot from 2024-06-26 13-55-07.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/c0c10fd0-eab1-4e63-8968-e8d50da268ebn%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/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/generic_linear_algebra.h>
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
// The program starts with the usual include files, all of which you should
// have seen before by now:
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/work_stream.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/fe/fe_dgq.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q.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>
// Then the usual placing of all content of this program into a namespace and
// the importation of the deal.II namespace into the one we will work in:
namespace Step26
{
using namespace dealii;
namespace EquationData
{
constexpr double eta = 1;
constexpr double kappa = 80.0; // nm²/ps
constexpr double beta = 10;
constexpr double density = 1;
template <int dim>
class TemperatureInitialValues : public Function<dim>
{
public:
TemperatureInitialValues()
: Function<dim>(1)
{}
virtual double value(const Point<dim> & /*p*/,
const unsigned int /*component*/ = 0) const override
{
return 1.0;
}
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override
{
for (unsigned int c = 0; c < this->n_components; ++c)
value(c) = TemperatureInitialValues<dim>::value(p, c);
}
};
template <int dim>
class TemperatureRightHandSide : public Function<dim>
{
public:
struct HeatSource
{
Point<dim> pos;
double E;
double radius;
double t;
};
TemperatureRightHandSide()
: Function<dim>()
, period(0.05)
{
Assert(dim == 2, ExcNotImplemented());
//read "heat_sources.txt" and parse points into a vector of vectors. Entries: x y z E radius t. Skip the first line.
std::ifstream file("heat_sources.txt");
AssertThrow(file.is_open(), ExcMessage("File not found"));
std::string line;
std::getline(file, line); //skip the first line
while (std::getline(file, line))
{
std::istringstream iss(line);
HeatSource heat_source;
iss >> heat_source.pos[0] >> heat_source.pos[1] >> heat_source.E >> heat_source.radius >> heat_source.t;
heat_sources.push_back(heat_source);
}
file.close();
}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override
{
(void)component;
AssertIndexRange(component, 1);
Assert(dim == 2, ExcNotImplemented());
double value = 0;
for (auto heat_source : active_sources)
{
const double distance = (p[0] - heat_source.pos[0]) * (p[0] - heat_source.pos[0]) +
(p[1] - heat_source.pos[1]) * (p[1] - heat_source.pos[1]);
if (distance < 2*heat_source.radius*heat_source.radius){
value += heat_source.E* std::exp(-distance / (2 * heat_source.radius * heat_source.radius));
}
}
return value;
}
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override
{
for (unsigned int c = 0; c < this->n_components; ++c)
value(c) = TemperatureRightHandSide<dim>::value(p, c);
}
bool is_source_active()
{
bool active = false;
active_sources.clear();
const double time = this->get_time();
for (const auto &heat_source : heat_sources)
{
const double t = heat_source.t;
if ((time >= t) && (time <= t + period)){
active_sources.push_back(heat_source);
active = true;
}
}
return active;
}
std::vector<HeatSource> active_sources;
private:
std::vector<HeatSource> heat_sources;
const double period;
};
}; // namespace EquationData
namespace Assembly
{
namespace Scratch
{
template <int dim>
struct TemperatureMatrix
{
TemperatureMatrix(const FiniteElement<dim> &fe,
const Mapping<dim> & mapping,
const Quadrature<dim> & temperature_quadrature);
TemperatureMatrix(const TemperatureMatrix &data);
FEValues<dim> fe_values;
std::vector<double> phi_T;
std::vector<Tensor<1, dim>> grad_phi_T;
};
template <int dim>
TemperatureMatrix<dim>::TemperatureMatrix(
const FiniteElement<dim> &fe,
const Mapping<dim> & mapping,
const Quadrature<dim> & temperature_quadrature)
: fe_values(mapping,
fe,
temperature_quadrature,
update_values | update_gradients |
update_JxW_values)
, phi_T(fe.n_dofs_per_cell())
, grad_phi_T(fe.n_dofs_per_cell())
{}
template <int dim>
TemperatureMatrix<dim>::TemperatureMatrix(
const TemperatureMatrix &scratch)
: fe_values(
scratch.fe_values.get_mapping(),
scratch.fe_values.get_fe(),
scratch.fe_values.get_quadrature(),
scratch.fe_values.get_update_flags())
, phi_T(scratch.phi_T)
, grad_phi_T(scratch.grad_phi_T)
{}
template <int dim>
struct TemperatureRHS
{
TemperatureRHS(const FiniteElement<dim> &fe,
const Mapping<dim> & mapping,
const Quadrature<dim> & quadrature);
TemperatureRHS(const TemperatureRHS &data);
FEValues<dim> fe_values;
std::vector<double> phi_T;
std::vector<Tensor<1, dim>> grad_phi_T;
std::vector<double> old_temperature_values;
std::vector<double> old_old_temperature_values;
std::vector<Tensor<1, dim>> old_temperature_grads;
std::vector<Tensor<1, dim>> old_old_temperature_grads;
std::vector<double> old_temperature_laplacians;
std::vector<double> old_old_temperature_laplacians;
};
template <int dim>
TemperatureRHS<dim>::TemperatureRHS(
const FiniteElement<dim> &fe,
const Mapping<dim> & mapping,
const Quadrature<dim> & quadrature)
: fe_values(mapping,
fe,
quadrature,
update_values | update_gradients |
update_hessians | update_quadrature_points |
update_JxW_values)
, phi_T(fe.n_dofs_per_cell())
, grad_phi_T(fe.n_dofs_per_cell())
,
old_temperature_values(quadrature.size())
, old_old_temperature_values(quadrature.size())
, old_temperature_grads(quadrature.size())
, old_old_temperature_grads(quadrature.size())
, old_temperature_laplacians(quadrature.size())
, old_old_temperature_laplacians(quadrature.size())
{}
template <int dim>
TemperatureRHS<dim>::TemperatureRHS(const TemperatureRHS &scratch)
: fe_values(
scratch.fe_values.get_mapping(),
scratch.fe_values.get_fe(),
scratch.fe_values.get_quadrature(),
scratch.fe_values.get_update_flags())
, phi_T(scratch.phi_T)
, grad_phi_T(scratch.grad_phi_T)
,
old_temperature_values(scratch.old_temperature_values)
, old_old_temperature_values(scratch.old_old_temperature_values)
, old_temperature_grads(scratch.old_temperature_grads)
, old_old_temperature_grads(scratch.old_old_temperature_grads)
, old_temperature_laplacians(scratch.old_temperature_laplacians)
, old_old_temperature_laplacians(scratch.old_old_temperature_laplacians)
{}
} // namespace Scratch
namespace CopyData
{
template <int dim>
struct TemperatureMatrix
{
TemperatureMatrix(const FiniteElement<dim> &fe);
FullMatrix<double> local_mass_matrix;
FullMatrix<double> local_stiffness_matrix;
std::vector<types::global_dof_index> local_dof_indices;
};
template <int dim>
TemperatureMatrix<dim>::TemperatureMatrix(
const FiniteElement<dim> &fe)
: local_mass_matrix(fe.n_dofs_per_cell(),
fe.n_dofs_per_cell())
, local_stiffness_matrix(fe.n_dofs_per_cell(),
fe.n_dofs_per_cell())
, local_dof_indices(fe.n_dofs_per_cell())
{}
template <int dim>
struct TemperatureRHS
{
TemperatureRHS(const FiniteElement<dim> &fe);
Vector<double> local_rhs;
std::vector<types::global_dof_index> local_dof_indices;
FullMatrix<double> matrix_for_bc;
};
template <int dim>
TemperatureRHS<dim>::TemperatureRHS(
const FiniteElement<dim> &fe)
: local_rhs(fe.n_dofs_per_cell())
, local_dof_indices(fe.n_dofs_per_cell())
, matrix_for_bc(fe.n_dofs_per_cell(),
fe.n_dofs_per_cell())
{}
} // namespace CopyData
} // namespace Assembly
template <int dim>
class HeatEquation
{
public:
HeatEquation();
void run();
private:
void setup_system();
void assemble_temperature_matrix();
void assemble_temperature_system();
void solve_time_step();
void output_results() const;
void refine_when_on(int max_grid_level);
void refine_mesh(const unsigned int max_grid_level);
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
const MappingQ<dim> mapping;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
AffineConstraints<double> constraints;
LA::MPI::Vector locally_relevant_solution;
LA::MPI::Vector locally_relevant_old_solution;
LA::MPI::Vector locally_relevant_old_old_solution;
LA::MPI::Vector system_rhs;
LA::MPI::SparseMatrix mass_matrix;
LA::MPI::SparseMatrix stiffness_matrix;
LA::MPI::SparseMatrix system_matrix;
std::shared_ptr<LA::MPI::PreconditionJacobi> T_preconditioner;
double time;
double time_step;
double old_time_step;
double t_off;
unsigned int timestep_number;
EquationData::TemperatureRightHandSide<dim> temperature_right_hand_side;
bool rebuild_temperature_matrices;
bool rebuild_temperature_preconditioner;
ConditionalOStream pcout;
TimerOutput computing_timer;
void setup_temperature_matrices(
const IndexSet &locally_owned_dofs,
const IndexSet &locally_relevant_dofs);
double get_cfl_number() const;
void local_assemble_temperature_matrix(
const typename DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::TemperatureMatrix<dim> & scratch,
Assembly::CopyData::TemperatureMatrix<dim> & data);
void copy_local_to_global_temperature_matrix(
const Assembly::CopyData::TemperatureMatrix<dim> &data);
void local_assemble_temperature_rhs(
const typename DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::TemperatureRHS<dim> & scratch,
Assembly::CopyData::TemperatureRHS<dim> & data);
void copy_local_to_global_temperature_rhs(
const Assembly::CopyData::TemperatureRHS<dim> &data);
class Postprocessor;
};
template <int dim>
HeatEquation<dim>::HeatEquation()
: mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator, typename Triangulation<dim>::MeshSmoothing(Triangulation<dim>::smoothing_on_refinement)) //| Triangulation<dim>::eliminate_refined_inner_islands))
, mapping(4)
, fe(1)
, dof_handler(triangulation)
, time_step(0)
, old_time_step(0)
, t_off(0)
, timestep_number(0)
, rebuild_temperature_matrices(true)
, rebuild_temperature_preconditioner(true)
, pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
, computing_timer(mpi_communicator, pcout, TimerOutput::never, TimerOutput::wall_times)
{}
template <int dim>
double HeatEquation<dim>::get_cfl_number() const
{
const QIterated<dim> quadrature_formula(QTrapezoid<1>(),
fe.degree);
const unsigned int n_q_points = quadrature_formula.size();
double max_local_cfl = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
for (unsigned int q = 0; q < n_q_points; ++q)
max_local_cfl =
std::max(max_local_cfl, cell->diameter());
//std::max(max_local_cfl, EquationData::kappa/(0.5*cell->diameter()*cell->diameter()));
}
return Utilities::MPI::max(max_local_cfl, MPI_COMM_WORLD);
}
template <int dim>
void HeatEquation<dim>::setup_temperature_matrices(
const IndexSet &locally_owned_dofs,
const IndexSet &locally_relevant_dofs)
{
T_preconditioner.reset();
mass_matrix.clear();
stiffness_matrix.clear();
system_matrix.clear();
DynamicSparsityPattern sp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern(dof_handler,
sp,
constraints,
false);
SparsityTools::distribute_sparsity_pattern(sp,
dof_handler.locally_owned_dofs(),
mpi_communicator,
locally_relevant_dofs);
sp.compress();
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
sp,
mpi_communicator);
//system_matrix.reinit(sp);
mass_matrix.reinit(system_matrix);
stiffness_matrix.reinit(system_matrix);
}
template <int dim>
void HeatEquation<dim>::setup_system()
{
TimerOutput::Scope timing_section(computing_timer, "Setup dof systems");
dof_handler.distribute_dofs(fe);
const types::global_dof_index n_T = dof_handler.n_dofs();
std::locale s = pcout.get_stream().getloc();
pcout.get_stream().imbue(std::locale(""));
pcout << "Number of active cells: " << triangulation.n_global_active_cells()
<< " (on " << triangulation.n_levels() << " levels)" << std::endl
<< "Number of degrees of freedom: " << n_T << ')' << std::endl
<< std::endl;
pcout.get_stream().imbue(s);
IndexSet locally_owned_dofs(n_T),
locally_relevant_dofs(n_T);
{
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
}
{
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler,
constraints);
VectorTools::interpolate_boundary_values(
dof_handler,
0,
EquationData::TemperatureInitialValues<dim>(),
constraints);
VectorTools::interpolate_boundary_values(
dof_handler,
1,
EquationData::TemperatureInitialValues<dim>(),
constraints);
constraints.close();
}
setup_temperature_matrices(locally_owned_dofs,
locally_relevant_dofs);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
locally_relevant_old_solution.reinit(locally_relevant_solution);
locally_relevant_old_old_solution.reinit(locally_relevant_solution);
rebuild_temperature_matrices = true;
rebuild_temperature_preconditioner = true;
}
template <int dim>
void HeatEquation<dim>::local_assemble_temperature_matrix(
const typename DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::TemperatureMatrix<dim> & scratch,
Assembly::CopyData::TemperatureMatrix<dim> & data)
{
const unsigned int dofs_per_cell =
scratch.fe_values.get_fe().n_dofs_per_cell();
const unsigned int n_q_points =
scratch.fe_values.n_quadrature_points;
scratch.fe_values.reinit(cell);
cell->get_dof_indices(data.local_dof_indices);
data.local_mass_matrix = 0;
data.local_stiffness_matrix = 0;
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
scratch.grad_phi_T[k] =
scratch.fe_values.shape_grad(k, q);
scratch.phi_T[k] = scratch.fe_values.shape_value(k, q);
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
data.local_mass_matrix(i, j) +=
(scratch.phi_T[i] * scratch.phi_T[j] *
scratch.fe_values.JxW(q));
data.local_stiffness_matrix(i, j) +=
(EquationData::kappa * scratch.grad_phi_T[i] *
scratch.grad_phi_T[j] * scratch.fe_values.JxW(q));
}
}
}
template <int dim>
void HeatEquation<dim>::copy_local_to_global_temperature_matrix(
const Assembly::CopyData::TemperatureMatrix<dim> &data)
{
constraints.distribute_local_to_global(data.local_mass_matrix,
data.local_dof_indices,
mass_matrix);
constraints.distribute_local_to_global(
data.local_stiffness_matrix,
data.local_dof_indices,
stiffness_matrix);
}
template <int dim>
void HeatEquation<dim>::assemble_temperature_matrix()
{
if (rebuild_temperature_matrices == false)
return;
TimerOutput::Scope timer_section(computing_timer,
" Assemble temperature matrices");
mass_matrix = 0;
stiffness_matrix = 0;
const QGauss<dim> quadrature_formula(fe.degree + 2);
using CellFilter =
FilteredIterator<typename DoFHandler<2>::active_cell_iterator>;
WorkStream::run(
CellFilter(IteratorFilters::LocallyOwnedCell(),
dof_handler.begin_active()),
CellFilter(IteratorFilters::LocallyOwnedCell(),
dof_handler.end()),
[this](const typename DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::TemperatureMatrix<dim> & scratch,
Assembly::CopyData::TemperatureMatrix<dim> & data) {
this->local_assemble_temperature_matrix(cell, scratch, data);
},
[this](const Assembly::CopyData::TemperatureMatrix<dim> &data) {
this->copy_local_to_global_temperature_matrix(data);
},
Assembly::Scratch::TemperatureMatrix<dim>(fe,
mapping,
quadrature_formula),
Assembly::CopyData::TemperatureMatrix<dim>(fe));
mass_matrix.compress(VectorOperation::add);
stiffness_matrix.compress(VectorOperation::add);
rebuild_temperature_matrices = false;
rebuild_temperature_preconditioner = true;
}
template <int dim>
void HeatEquation<dim>::local_assemble_temperature_rhs(
const typename DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::TemperatureRHS<dim> & scratch,
Assembly::CopyData::TemperatureRHS<dim> & data)
{
const bool use_bdf2_scheme = (timestep_number != 0);
const unsigned int dofs_per_cell =
scratch.fe_values.get_fe().n_dofs_per_cell();
const unsigned int n_q_points =
scratch.fe_values.n_quadrature_points;
data.local_rhs = 0;
cell->get_dof_indices(data.local_dof_indices);
scratch.fe_values.reinit(cell);
scratch.fe_values.get_function_values(
locally_relevant_old_solution, scratch.old_temperature_values);
scratch.fe_values.get_function_values(
locally_relevant_old_old_solution, scratch.old_old_temperature_values);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
scratch.phi_T[k] = scratch.fe_values.shape_value(k, q);
scratch.grad_phi_T[k] =
scratch.fe_values.shape_grad(k, q);
}
const double T_term_for_rhs =
(use_bdf2_scheme ?
(scratch.old_temperature_values[q] *
(1 + time_step / old_time_step) -
scratch.old_old_temperature_values[q] * (time_step * time_step) /
(old_time_step * (time_step + old_time_step))) :
scratch.old_temperature_values[q]);
const double gamma = temperature_right_hand_side.value(scratch.fe_values.quadrature_point(q), 0);
if (gamma>0.){
t_off = 0;
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
data.local_rhs(i) += (T_term_for_rhs + time_step * gamma ) * scratch.fe_values.JxW(q) * scratch.phi_T[i];
//pcout << "asdASDEHHD"<< constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]) << std::endl;
if (constraints.is_inhomogeneously_constrained(
data.local_dof_indices[i]))
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
data.matrix_for_bc(j, i) +=0 ;
/*
(scratch.phi_T[i] * scratch.phi_T[j] *
(use_bdf2_scheme ? ((2 * time_step + old_time_step) /
(time_step + old_time_step)) :
1.) ) *
scratch.fe_values.JxW(q);
*/
}
}
}
}
template <int dim>
void HeatEquation<dim>::copy_local_to_global_temperature_rhs(
const Assembly::CopyData::TemperatureRHS<dim> &data)
{
constraints.distribute_local_to_global(data.local_rhs,
data.local_dof_indices,
system_rhs,
data.matrix_for_bc);
}
template <int dim>
void HeatEquation<dim>::assemble_temperature_system()
{
const bool use_bdf2_scheme = (timestep_number != 0);
if (use_bdf2_scheme == true)
{
system_matrix.copy_from(mass_matrix);
system_matrix *=
(2 * time_step + old_time_step) / (time_step + old_time_step);
system_matrix.add(time_step, stiffness_matrix);
}
else
{
system_matrix.copy_from(mass_matrix);
system_matrix.add(time_step, stiffness_matrix);
}
if (rebuild_temperature_preconditioner == true)
{
T_preconditioner = std::make_shared<LA::MPI::PreconditionJacobi>();
T_preconditioner->initialize(system_matrix);
rebuild_temperature_preconditioner = false;
}
// The next part is computing the right hand side vectors. To do so, we
// first compute the average temperature $T_m$ that we use for evaluating
// the artificial viscosity stabilization through the residual $E(T) =
// (T-T_m)^2$. We do this by defining the midpoint between maximum and
// minimum temperature as average temperature in the definition of the
// entropy viscosity. An alternative would be to use the integral average,
// but the results are not very sensitive to this choice. The rest then
// only requires calling WorkStream::run again, binding the arguments to
// the <code>local_assemble_system_rhs</code> function that are the
// same in every call to the correct values:
system_rhs = 0;
const QGauss<dim> quadrature_formula(fe.degree + 2);
using CellFilter =
FilteredIterator<typename DoFHandler<2>::active_cell_iterator>;
auto worker =
[this](
const typename DoFHandler<dim>::active_cell_iterator &cell,
Assembly::Scratch::TemperatureRHS<dim> & scratch,
Assembly::CopyData::TemperatureRHS<dim> & data) {
this->local_assemble_temperature_rhs(cell,
scratch,
data);
};
auto copier = [this](const Assembly::CopyData::TemperatureRHS<dim> &data) {
this->copy_local_to_global_temperature_rhs(data);
};
WorkStream::run(CellFilter(IteratorFilters::LocallyOwnedCell(),
dof_handler.begin_active()),
CellFilter(IteratorFilters::LocallyOwnedCell(),
dof_handler.end()),
worker,
copier,
Assembly::Scratch::TemperatureRHS<dim>(
fe, mapping, quadrature_formula),
Assembly::CopyData::TemperatureRHS<dim>(fe));
system_rhs.compress(VectorOperation::add);
}
template <int dim>
void HeatEquation<dim>::solve_time_step()
{
{
TimerOutput::Scope timer_section(computing_timer,
" Assemble temperature rhs");
old_time_step = time_step;
const double scaling = (dim == 3 ? 0.25 : 1.0);
time_step = 1/(200.);//*0.353553/get_cfl_number(); //scaling / (2.1 * dim * std::sqrt(1. * dim)) /(fe.degree * get_cfl_number());
locally_relevant_solution = locally_relevant_old_solution;
assemble_temperature_system();
}
{
TimerOutput::Scope timer_section(computing_timer,
" Solve temperature system");
SolverControl solver_control(system_matrix.m(),
1e-12 * system_rhs.l2_norm());
SolverCG<LA::MPI::Vector> cg(solver_control);
LA::MPI::Vector distributed_temperature_solution(system_rhs);
distributed_temperature_solution = locally_relevant_solution;
cg.solve(system_matrix,
distributed_temperature_solution,
system_rhs,
*T_preconditioner);
constraints.distribute(distributed_temperature_solution);
locally_relevant_solution = distributed_temperature_solution;
pcout << " " << solver_control.last_step()
<< " CG iterations for temperature" << std::endl;
double temperature[2] = {std::numeric_limits<double>::max(),
-std::numeric_limits<double>::max()};
double global_temperature[2];
for (unsigned int i =
distributed_temperature_solution.local_range().first;
i < distributed_temperature_solution.local_range().second;
++i)
{
temperature[0] =
std::min<double>(temperature[0],
distributed_temperature_solution(i));
temperature[1] =
std::max<double>(temperature[1],
distributed_temperature_solution(i));
}
temperature[0] *= -1.0;
Utilities::MPI::max(temperature, MPI_COMM_WORLD, global_temperature);
global_temperature[0] *= -1.0;
pcout << " Temperature range: " << global_temperature[0] << ' '
<< global_temperature[1] << std::endl;
}
}
template <int dim>
void HeatEquation<dim>::refine_when_on( int max_grid_level)
{
parallel::distributed::SolutionTransfer<dim, LA::MPI::Vector> temperature_trans(dof_handler);
int n = 0;
while (n< max_grid_level){
for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(); cell != triangulation.end(); ++cell){
cell->clear_refine_flag();
if(cell->level() >= max_grid_level){
continue;
}
for (auto source : temperature_right_hand_side.active_sources){
if ((cell->center()-source.pos)*(cell->center()-source.pos) < 4*source.radius*source.radius || cell->point_inside(source.pos)){
cell->set_refine_flag();
}
}
}
std::vector<const LA::MPI::Vector *> x_temperature(2);
x_temperature[0] = &locally_relevant_solution;
x_temperature[1] = &locally_relevant_old_solution;
triangulation.prepare_coarsening_and_refinement();
temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
triangulation.execute_coarsening_and_refinement();
setup_system();
{
{
LA::MPI::Vector distributed_temp1(system_rhs);
LA::MPI::Vector distributed_temp2(system_rhs);
std::vector<LA::MPI::Vector *> tmp(2);
tmp[0] = &(distributed_temp1);
tmp[1] = &(distributed_temp2);
temperature_trans.interpolate(tmp);
constraints.distribute(distributed_temp1);
constraints.distribute(distributed_temp2);
locally_relevant_solution = distributed_temp1;
locally_relevant_old_solution = distributed_temp2;
}
}
n+=1;
}
}
template <int dim>
void HeatEquation<dim>::refine_mesh(unsigned int max_grid_level)
{
parallel::distributed::SolutionTransfer<dim, LA::MPI::Vector> temperature_trans(dof_handler);
{
TimerOutput::Scope timer_section(computing_timer,
"Refine mesh structure, part 1");
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> *>(),
locally_relevant_solution,
estimated_error_per_cell,
ComponentMask(),
nullptr,
0,
triangulation.locally_owned_subdomain());
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
triangulation, estimated_error_per_cell, 0.6, 0.4);
if (t_off > 0.25 && t_off <= 0.5){
max_grid_level -= 1;
}
else if(t_off > 0.5 && t_off <= 0.75){
max_grid_level -= 2;
}
else if(t_off > 0.75 && t_off <= 1.){
max_grid_level -= 3;
}
else if(t_off > 1.){
max_grid_level -= 4;
}
if (triangulation.n_levels() > max_grid_level)
for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(max_grid_level);
cell != triangulation.end();
++cell)
cell->clear_refine_flag();
std::vector<const LA::MPI::Vector *> x_temperature(2);
x_temperature[0] = &locally_relevant_solution;
x_temperature[1] = &locally_relevant_old_solution;
triangulation.prepare_coarsening_and_refinement();
temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
triangulation.execute_coarsening_and_refinement();
}
setup_system();
{
TimerOutput::Scope timer_section(computing_timer,
"Refine mesh structure, part 2");
{
LA::MPI::Vector distributed_temp1(system_rhs);
LA::MPI::Vector distributed_temp2(system_rhs);
std::vector<LA::MPI::Vector *> tmp(2);
tmp[0] = &(distributed_temp1);
tmp[1] = &(distributed_temp2);
temperature_trans.interpolate(tmp);
constraints.distribute(distributed_temp1);
constraints.distribute(distributed_temp2);
locally_relevant_solution = distributed_temp1;
locally_relevant_old_solution = distributed_temp2;
}
}
}
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(locally_relevant_solution, "U");
Vector<float> subdomain(triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches();
//data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number));
const std::string filename =
"solution-";
//"solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
static int out_index = 0;
//data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
data_out.write_vtu_with_pvtu_record("./", filename, out_index, mpi_communicator, 2, 1);
out_index++;
}
template <int dim>
void HeatEquation<dim>::run()
{
const unsigned int initial_global_refinement = (dim == 2 ? 4 : 3);
const unsigned int initial_adaptive_refinement = (dim == 2 ? 4 : 2);
const unsigned int adaptive_refinement_interval = 10;
bool generate_graphical_output = true;
int graphical_output_interval = 1;
double end_time = 1.0;
bool source_on = false;
//GridGenerator::hyper_L(triangulation);
GridGenerator::hyper_cube( triangulation, 0, 300, false);
triangulation.refine_global(initial_global_refinement);
setup_system();
unsigned int pre_refinement_step = 0;
start_time_iteration:
{
LA::MPI::Vector solution(dof_handler.locally_owned_dofs(), mpi_communicator);
VectorTools::project(dof_handler,
constraints,
QGauss<dim>(fe.degree + 2),
EquationData::TemperatureInitialValues<dim>(),
solution);
locally_relevant_solution = solution;
locally_relevant_old_solution = solution;
locally_relevant_old_old_solution = solution;
}
timestep_number = 0;
time_step = old_time_step = 0;
double time = 0;
temperature_right_hand_side.set_time(time);
do
{
pcout << "Timestep " << timestep_number
<< ": t=" << time << std::endl;
temperature_right_hand_side.set_time(time);
bool is_source_active = temperature_right_hand_side.is_source_active();
assemble_temperature_matrix();
solve_time_step();
pcout << std::endl;
if ((timestep_number == 0) &&
(pre_refinement_step < initial_adaptive_refinement))
{
if (is_source_active && source_on == false)
{
source_on = true;
refine_when_on(initial_global_refinement + initial_adaptive_refinement);
}
else if (!is_source_active && source_on == true){ source_on = false; }
refine_when_on(initial_global_refinement + initial_adaptive_refinement);
refine_mesh(initial_global_refinement + initial_adaptive_refinement);
++pre_refinement_step;
goto start_time_iteration;
}
if (is_source_active && source_on == false)
{
source_on = true;
refine_when_on(initial_global_refinement + initial_adaptive_refinement);
}
else if (!is_source_active && source_on == true){source_on = false;}
if ((timestep_number > 0) && (timestep_number % adaptive_refinement_interval == 0))
{
refine_mesh(initial_global_refinement +
initial_adaptive_refinement);
}
if ((generate_graphical_output == true) &&
(timestep_number % graphical_output_interval == 0))
output_results();
if (time > end_time )
break;
locally_relevant_old_old_solution = locally_relevant_old_solution;
locally_relevant_old_solution = locally_relevant_solution;
if (old_time_step > 0)
{
{
LA::MPI::Vector distr_solution(system_rhs);
distr_solution = locally_relevant_solution;
LA::MPI::Vector distr_old_solution(system_rhs);
distr_old_solution = locally_relevant_old_old_solution;
distr_solution.sadd(1. + time_step / old_time_step,
-time_step / old_time_step,
distr_old_solution);
locally_relevant_solution = distr_solution;
}
}
if ((timestep_number > 0) && (timestep_number % 100 == 0))
computing_timer.print_summary();
time += time_step;
//pcout << "t_off = " << t_off << std::endl;
t_off += time_step;
++timestep_number;
//if (timestep_number == 1) break; // Flag for debugging
}
while (true);
// If we are generating graphical output, do so also for the last time
// step unless we had just done so before we left the do-while loop
if ((generate_graphical_output == true) &&
!((timestep_number - 1) % graphical_output_interval == 0))
output_results();
}
} // namespace Step26
int main(int argc, char *argv[])
{
try
{ using namespace dealii;
using namespace Step26;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
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;
}
