Hello,
I've been working to modify the Step 26 tutorial to include a radiative
boundary condition on one edge and am running into an issue combining this
with adaptive mesh refinement. Although the boundary condition and adaptive
refinement both seem to be working separately, I get large artifacts at
the radiative boundary after each timestep where I refine the mesh.
During my attempts to debug I've found several other strange symptoms:
1) This problem doesn't occur after I initially refine the mesh (presumably
since I set it to a predefined initial temperature?)
2) The artifacts don't appear if I run the code without ever calling
refine_mesh() (near-identical to the function by the same name in the
tutorial code), but they *do* appear if I call refine_mesh() but set the
refine/coarsen fractions to 0 so no refinement actually occurs.
3) If I set a fixed "local" temperature instead of using
get_function_values to estimate a local temperature at each quadrature
point the problem goes away. This is the most mysterious to me, since it
seems to suggest I might be accessing the local temperatures incorrectly
(and the temperature values at the surface before/after the refinement seem
consistent with this) but I haven't been able to find a bug there yet.
Does any of this suggest something obvious I might be missing? I'm still
pretty new to deal.II so I wouldn't be surprised if the bug was something
very simple.
In case it's helpful, I've included a (still kind of long) minimal version
of the code I'm using as well as images showing examples of a solution
before and after the refinement step (note that the refine/coarsen
fractions are zero in these images so the mesh is not actually being
refined).
Best,
Sarah
--
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/de0e5cf9-650a-4082-a70b-ea04e5cae856n%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2019 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/exceptions.h>
#include <deal.II/base/function.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.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_update_flags.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/data_out_faces.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <cstring>
#include <exception>
#include <iostream>
#include <filesystem>
#include <map>
#include <string>
#include <vector>
#include <cmath>
#include "../support_code/config_in.cc"
#include "../support_code/timer.h"
#include "heat_transfer.cc"
#include "rheology.cc"
#include "heat_right_hand_side.cc"
#include "initial_temperature.cc"
#include "heat_boundary_values.cc"
#include <armadillo>
using namespace dealii;
using namespace arma;
namespace fs = std::filesystem;
template <int dim>
class HeatEquation
{
public:
HeatEquation(config_in&);
void clear_output_directory();
void run();
config_in& cfg;
private:
void load_initial_temperature();
void heat_setup_system();
void heat_compute_mass_and_laplace_matrices();
void heat_setup_crank_nicolson();
void heat_solve_time_step();
void do_heat_step();
void refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level);
void graphical_output_results() const;
void textual_output_results() const;
void textual_output_results_blocks() const;
void set_initial_temperature();
void set_boundary_indicators();
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> heat_mass_matrix;
SparseMatrix<double> heat_bc_matrix;
SparseMatrix<double> heat_laplace_matrix;
SparseMatrix<double> heat_system_matrix;
Vector<double> heat_tmp;
Vector<double> heat_bc_tmp;
Vector<double> heat_forcing_terms;
Vector<double> heat_solution;
Vector<double> old_heat_solution;
Vector<double> heat_system_rhs;
Vector<double> heat_bc_rhs;
LinearAlgebra::Vector<double> heat_solution_la;
HeatBoundaryValuesRight<dim> heat_boundary_values_function_right;
HeatBoundaryValuesBottom<dim> heat_boundary_values_function_bottom;
double time;
double time_step;
unsigned int timestep_number;
const double theta;
vec x_vec;
vec z_vec;
mat initial_temperature_mat;
double T_eq;
};
template <int dim>
HeatEquation<dim>::HeatEquation(config_in& cfgi)
: fe(1)
, dof_handler(triangulation)
, cfg(cfgi)
, time(cfgi.present_time)
, time_step(cfgi.time_step)
, timestep_number(0)
, theta(0.5)
{};
template <int dim>
void HeatEquation<dim>::load_initial_temperature()
{
x_vec.load(cfg.x_file, raw_ascii);
z_vec.load(cfg.z_file, raw_ascii);
initial_temperature_mat.load(cfg.temp_file, raw_ascii);
T_eq = cfg.T_surf;
}
template <int dim>
void HeatEquation<dim>::heat_setup_system()
{
dof_handler.distribute_dofs(fe);
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);
heat_mass_matrix.reinit(sparsity_pattern);
heat_laplace_matrix.reinit(sparsity_pattern);
heat_bc_matrix.reinit(sparsity_pattern);
heat_system_matrix.reinit(sparsity_pattern);
heat_solution.reinit(dof_handler.n_dofs());
old_heat_solution.reinit(dof_handler.n_dofs());
heat_bc_rhs.reinit(dof_handler.n_dofs());
heat_system_rhs.reinit(dof_handler.n_dofs());
if (timestep_number == 0)
{
set_initial_temperature();
std::cout << "Set initial temperature" << std::endl;
}
}
template <int dim>
void HeatEquation<dim>::heat_compute_mass_and_laplace_matrices()
{
double T_local;
const QGauss<dim> quadrature_formula(fe.degree + 1);
const QGauss<dim-1> top_quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_top_values(fe,
top_quadrature_formula,
update_values |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_top_q_points = top_quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<double> temperature_at_q_points(n_q_points);
std::vector<double> temperature_at_top_q_points(n_top_q_points);
// make mass and laplace matrices
MatrixCreator::create_mass_matrix(dof_handler,
QGauss<dim>(fe.degree + 1),
heat_mass_matrix);
heat_mass_matrix *= 1000*3000;
MatrixCreator::create_laplace_matrix(dof_handler,
QGauss<dim>(fe.degree + 1),
heat_laplace_matrix);
heat_laplace_matrix *= 1e-6;
// make matrix + rhs contributions from radiative boundary
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
fe_values.get_function_values(heat_solution,temperature_at_q_points);
for (unsigned int f = 0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
if (cell->face(f)->at_boundary() && (cell->face(f)->boundary_id() == 3))
{
fe_top_values.reinit(cell, f);
fe_top_values.get_function_values(heat_solution,temperature_at_top_q_points); //might be the problem?
for (unsigned int q_index = 0; q_index < n_top_q_points; ++q_index)
{
T_local = temperature_at_top_q_points[q_index];
std::cout << T_local << std::endl;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) += (5e-8 * // sigma
pow(T_eq,4) * //T^{n-1}^4
fe_top_values.shape_value(i,q_index) * //
fe_top_values.JxW(q_index));
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
(5e-8 * pow(T_local,3) * // a(x_q)
fe_top_values.shape_value(i, q_index) * // phi_i(x_q)
fe_top_values.shape_value(j, q_index) * // phi_j(x_q)
fe_top_values.JxW(q_index)); // dx
}
}
}
}
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
heat_bc_matrix.add(local_dof_indices[i], local_dof_indices[j], cell_matrix(i, j));
heat_bc_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
}
template <int dim>
void HeatEquation<dim>::heat_setup_crank_nicolson()
{
heat_tmp.reinit(heat_solution.size());
heat_bc_tmp.reinit(heat_solution.size());
heat_forcing_terms.reinit(heat_solution.size());
heat_mass_matrix.vmult(heat_system_rhs, old_heat_solution);
heat_laplace_matrix.vmult(heat_tmp, old_heat_solution);
heat_bc_matrix.vmult(heat_bc_tmp, old_heat_solution);
heat_system_rhs.add(-(1 - theta) * time_step, heat_tmp);
heat_system_rhs.add((1 - 4*theta) * time_step , heat_bc_tmp);
heat_system_rhs.add(-time_step, heat_bc_rhs);
heat_system_matrix.copy_from(heat_mass_matrix);
heat_system_matrix.add(theta * time_step, heat_laplace_matrix);
heat_system_matrix.add(-4 * theta * time_step, heat_bc_matrix);
constraints.condense(heat_system_matrix, heat_system_rhs);
{
std::map<types::global_dof_index, double> heat_boundary_values_right;
std::map<types::global_dof_index, double> heat_boundary_values_bottom;
VectorTools::interpolate_boundary_values(dof_handler,
1,
heat_boundary_values_function_bottom,
heat_boundary_values_bottom);
VectorTools::interpolate_boundary_values(dof_handler,
2,
heat_boundary_values_function_right,
heat_boundary_values_right);
MatrixTools::apply_boundary_values(heat_boundary_values_bottom,
heat_system_matrix,
heat_solution,
heat_system_rhs);
MatrixTools::apply_boundary_values(heat_boundary_values_right,
heat_system_matrix,
heat_solution,
heat_system_rhs);
}
}
template <int dim>
void HeatEquation<dim>::heat_solve_time_step()
{
SolverControl solver_control(cfg.heat_iteration_coefficient, cfg.heat_tolerance_coefficient * heat_system_rhs.l2_norm());
SolverCG<> cg(solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(heat_system_matrix, 1.0);
cg.solve(heat_system_matrix, heat_solution, heat_system_rhs, preconditioner);
constraints.distribute(heat_solution);
std::cout << " " << solver_control.last_step() << " CG iterations."
<< std::endl;
}
template <int dim>
void HeatEquation<dim>::do_heat_step()
{
heat_compute_mass_and_laplace_matrices();
heat_setup_crank_nicolson();
heat_solve_time_step();
graphical_output_results();
}
template <int dim>
void HeatEquation<dim>::graphical_output_results() const
{
if (timestep_number % 1 == 0) {
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(heat_solution, "U");
data_out.add_data_vector(old_heat_solution, "U_old");
data_out.build_patches();
const std::string filename = cfg.output_folder +
"/solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
std::ofstream output(filename);
data_out.write_vtk(output);
}
}
template <int dim>
void HeatEquation<dim>::set_initial_temperature()
{
InitialTemperature<dim> t_init;
t_init.set_initial_temperature_field(x_vec,z_vec,initial_temperature_mat);
VectorTools::interpolate(dof_handler,
t_init,
old_heat_solution);
heat_solution = old_heat_solution;
}
template <int dim>
void HeatEquation<dim>::set_boundary_indicators()
{
double zero_tolerance = 1.0e-4;
double distance_left;
double distance_right;
double distance_bottom;
double distance_top;
for (const auto &cell : dof_handler.active_cell_iterators()) // loop over all cells
{
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) // loop over all faces
{
if (cell->face(f)->at_boundary())
{
distance_left = fabs(cell->face(f)->center()[0] - cfg.x_left);
distance_right = fabs(cell->face(f)->center()[0] - cfg.x_right);
distance_bottom = fabs(cell->face(f)->center()[1] - cfg.z_bottom);
distance_top = fabs(cell->face(f)->center()[1]);
if (distance_left < zero_tolerance)
{
cell->face(f)->set_all_boundary_ids(0);
}
else if (distance_bottom < zero_tolerance)
{
cell->face(f)->set_all_boundary_ids(1);
}
else if (distance_right < zero_tolerance)
{
cell->face(f)->set_all_boundary_ids(2);
}
else if (distance_top < zero_tolerance)
{
cell->face(f)->set_all_boundary_ids(3);
}
}
}
}
}
template <int dim>
void HeatEquation<dim>::refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
Vector<float> gradient_indicator(triangulation.n_active_cells());
DerivativeApproximation::approximate_gradient(dof_handler,
heat_solution,
gradient_indicator);
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
gradient_indicator,
0.1,
0.15);
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 = heat_solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
triangulation.execute_coarsening_and_refinement();
heat_setup_system();
solution_trans.interpolate(previous_solution, heat_solution);
constraints.distribute(heat_solution);
}
template <int dim>
void HeatEquation<dim>::run()
{
clear_output_directory();
GridIn<dim> grid_in;
grid_in.attach_triangulation(triangulation);
std::ifstream mesh_stream(cfg.mesh_filename,
std::ifstream::in);
grid_in.read_ucd(mesh_stream);
set_boundary_indicators();
load_initial_temperature();
heat_setup_system();
unsigned int pre_refinement_step = 0;
start_time_iteration:
time = 0.0;
timestep_number = 0;
heat_tmp.reinit(heat_solution.size());
heat_bc_tmp.reinit(heat_solution.size());
heat_forcing_terms.reinit(heat_solution.size());
heat_boundary_values_function_right.set_temperature_field(x_vec,z_vec,initial_temperature_mat);
heat_boundary_values_function_bottom.set_temperature_field(x_vec,z_vec,initial_temperature_mat);
graphical_output_results();
set_boundary_indicators();
heat_setup_system();
while (time <= cfg.final_time)
{
time += time_step;
++timestep_number;
std::cout << "Time step " << timestep_number << std::endl;
std::cout << " at t= " << time/SECSINYEAR/1e6 << " My" << std::endl;
do_heat_step();
if ((timestep_number == 1) && (pre_refinement_step < cfg.adaptive_refinement))
{
refine_mesh(cfg.global_refinement,
cfg.global_refinement +
cfg.adaptive_refinement);
++pre_refinement_step;
heat_tmp.reinit(heat_solution.size());
heat_bc_tmp.reinit(heat_solution.size());
heat_forcing_terms.reinit(heat_solution.size());
std::cout << std::endl;
goto start_time_iteration;
}
else if ((timestep_number > 0) && (timestep_number % 5 == 0))
{
refine_mesh(cfg.global_refinement,
cfg.global_refinement +
cfg.adaptive_refinement);
heat_tmp.reinit(heat_solution.size());
heat_bc_tmp.reinit(heat_solution.size());
heat_forcing_terms.reinit(heat_solution.size());
}
old_heat_solution = heat_solution;
}
}
int main(int argc, char* argv[])
{
try
{
using namespace dealii;
char* config_filename = new char[120];
if (argc == 1) // if no input parameters (as if launched from eclipse)
{
std::strcpy(config_filename,"config.cfg");
}
config_in cfg(config_filename);
HeatEquation<2> heat_equation_solver(cfg);
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;
}
