This question might be related to earlier of my questions, but finally I could not find an explanation for the behaviour yet, thus this question.
I intended to calculate the time-dependent heat equation with constant factors, using the newton method (such that time-dependent factors can be implemented later) and automatic differentiation. Now I added an offset to the initial values and the border values, such that both are lifted equally. Nevertheless I found out that if I either reduce the time step length or increase the offset value, the residual value I am using for controlling the solution progress initially decreases fast, but then slows down at a certain value depending on the offset value and the time step size. Below a certain time step size or above a certain offset value the value never gets smaller, and the best calculated step length goes towards 0 (compared to close to 1 before). I do not understand that behaviour. Based on initial tests my code should work correct, nevertheless it quits for certain values. What could be the reason for that? I attached the code to the question, it should be self-contained. Changing of the parameters can either be done using the OFFSET-parameter or the time_step parameter in the parameter.prm-file Thanks! -- 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 dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/data_out_base.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/solver_qmrs.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_sparsity_pattern.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/vector_tools.templates.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <Sacado.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/filesystem.hpp>
#include <boost/mpi.hpp>
#include <boost/mpi/collectives.hpp>
#include <boost/serialization/serialization.hpp>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <experimental/type_traits>
#include <map>
#include <utility>
#include <armadillo>
#define INNER_CIRCLE_RADIUS 0.1
#define PARAMETER_VERSION 2
#define USE_BROADCAST
#define STATIC_HEAT_SOURCE_VALUE_OUTER 0
#define STATIC_HEAT_SOURCE_VALUE_INNER 1e14
#define OFFSET 1e2
template < typename T >
using toDim_t = decltype(std::declval<T>().norm());
template < typename T >
using has_toDim = std::experimental::is_detected< toDim_t, T >;
template <bool, typename T = void>
struct enable_if
{};
template <typename T>
struct enable_if<true, T> {
typedef T type;
};
namespace Step15
{
using namespace dealii;
void print_status_update(const ConditionalOStream &pcout, const std::string input, const bool do_it = false)
{
#ifndef NO_COMMENTS
if(do_it)
{
pcout << input << '\n';
//getchar();
}
#endif
}
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues() : Function<dim>(1 * dim) {}
virtual double value(const Point<dim> &p, const unsigned int component) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &values) const;
};
template <int dim>
double InitialValues<dim>::value(const Point<dim> &p, const unsigned int component) const
{
(void) p;
(void) component;
return OFFSET;
}
template <int dim>
void InitialValues<dim>::vector_value(const Point<dim> &p, Vector<double> &values) const
{
for(size_t i = 0; i < values.size(); ++i)
values[i] = value(p, i);
}
template <int dim>
class MinimalSurfaceProblem
{
public:
MinimalSurfaceProblem (const char * file_name);
~MinimalSurfaceProblem ();
void run ();
private:
void setup_system (const bool initial_step);
void setup_matrices(const IndexSet &partitioner, const IndexSet &relevant_partitioner);
void assemble_system (void);
void solve (void);
void refine_mesh ();
void set_boundary_values ();
double compute_residual (const double alpha);
double determine_step_length ();
double recalculate_step_length(void);
void output_results(const int cycle) const;
void print_usage_message(void);
void declare_parameters(void);
void parse_parameters(void);
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
const size_t n_components;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
ConstraintMatrix hanging_node_constraints;
ConstraintMatrix newton_constraints;
ConstraintMatrix boundary_constraints;
LinearAlgebraTrilinos::MPI::SparseMatrix system_matrix;
LinearAlgebraTrilinos::MPI::Vector present_solution;
LinearAlgebraTrilinos::MPI::Vector old_solution;
LinearAlgebraTrilinos::MPI::Vector newton_update;
LinearAlgebraTrilinos::MPI::Vector system_rhs;
LinearAlgebraTrilinos::MPI::Vector residual;
LinearAlgebraTrilinos::MPI::Vector evaluation_point;
LinearAlgebraTrilinos::MPI::Vector extension_vector;
std::string log_file_name;
std::ostream & out_stream;
std::ofstream file_out_stream;
std::string data_path;
ConditionalOStream pcout;
ParameterHandler prm;
size_t max_grid_level;
size_t min_grid_level;
size_t max_inner_iterations;
size_t max_same_res_values;
double start_time;
double max_time;
double time_step;
double theta;
double alpha_val;
arma::mat intensity_matrix, old_intensity_matrix;
double pulse_energy;
int pulse_calculation_points;
double pulse_duration;
double beam_width;
double focal_length;
double focal_length_inside;
double pulse_wavelength;
std::string parameters_file_name;
std::vector<bool> boundary_conditions; //Order of processing: Up, down, left, right. True: Dirichlet, False: Neumann
const size_t Electron_Temperature = 0;
};
// @sect3{Boundary condition}
template <int dim>
class BoundaryValuesCarriers : public Function<dim>
{
public:
BoundaryValuesCarriers (const size_t n_components) : Function<dim>(1 * dim), n_components(n_components) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &value) const;
private:
const size_t n_components;
};
template <int dim>
double BoundaryValuesCarriers<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
(void) component;
(void) p;
return 0;
//return std::sin(2 * numbers::PI * (p[0]+p[1]));
}
template <int dim>
void BoundaryValuesCarriers<dim>::vector_value(const Point<dim> &p, Vector<double> &value) const
{
for(size_t i = 0; i < value.size(); ++i)
value[i] = BoundaryValuesCarriers<dim>::value(p, i);
}
template <int dim>
class BoundaryValuesTemperatures : public Function<dim>
{
public:
BoundaryValuesTemperatures (const size_t n_components) : Function<dim>(dim), n_components(n_components) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &value) const;
private:
const size_t n_components;
};
template <int dim>
double BoundaryValuesTemperatures<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
(void) component;
return p[0] * p[0] + OFFSET;
}
template <int dim>
void BoundaryValuesTemperatures<dim>::vector_value(const Point<dim> &p, Vector<double> &value) const
{
for(size_t i = 0; i < value.size(); ++i)
value[i] = BoundaryValuesTemperatures<dim>::value(p, i);
}
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues (const size_t n_components) : Function<dim>(dim), n_components(n_components) {}
virtual double value (const Point<dim> &p,
const size_t component = 0) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &value) const;
private:
const size_t n_components;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> &p, const size_t component) const
{
return BoundaryValuesTemperatures<dim>::value(p, component);
}
template <int dim>
void BoundaryValues<dim>::vector_value(const Point<dim> &p, Vector<double> &value) const
{
return BoundaryValuesTemperatures<dim>::vector_value(p, value);
}
template <int dim, typename R, typename T, typename U, typename V>
T calculate_right_side(const T &val_TE, const T &val_TE_old, const U &grad_TE, const U &grad_TE_old,
const Tensor<2, dim> &TE_gradient, const V &JxW_val,
const double theta)
{
(void) val_TE;
(void) val_TE_old;
R return_value = 0;
for(size_t d = 0; d < dim; ++d)
{
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
return_value += (theta * grad_TE[d][local_dim]
+ (1 - theta) * grad_TE_old[d][local_dim])
* TE_gradient[d][local_dim]
* JxW_val;
}
}
return return_value;
}
template <int dim>
MinimalSurfaceProblem<dim>::MinimalSurfaceProblem (const char *file_name)
:
mpi_communicator(MPI_COMM_WORLD),
triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(Triangulation<dim>::smoothing_on_refinement | Triangulation<dim>::smoothing_on_coarsening)),
n_components(1),
dof_handler (triangulation),
fe (FE_Q<dim>(2), n_components * dim),
out_stream(std::cout),
pcout(out_stream, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
alpha_val(0.1)
{
parameters_file_name = std::string(file_name);
this->declare_parameters();
prm.parse_input(file_name);
this->parse_parameters();
}
template <int dim>
MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem ()
{
file_out_stream.close();
dof_handler.clear ();
}
template <int dim>
void MinimalSurfaceProblem<dim>::print_usage_message()
{
static const char *message
=
"\n"
"TTM-calculator for SI.\n"
"\n"
"Usage:\n"
" ./main [-p parameter_file]\n"
"\n"
"The parameter file has the following format and allows the following\n"
"values (you can cut and paste this and use it for your own parameter\n"
"file):\n"
"\n";
std::cout << message;
prm.print_parameters (std::cout, ParameterHandler::Text);
}
template <int dim>
void MinimalSurfaceProblem<dim>::declare_parameters()
{
prm.declare_entry("Version", "0", Patterns::Integer(), "Sets the version of the parameter file, allows determination if the parameter file is new enough");
prm.enter_subsection("File processing data");
{
prm.declare_entry("Use stdcout", "true", Patterns::Bool(), "If set to true, stdcout will be used, else a file will be used");
prm.declare_entry("Log file name", "log_file.txt", Patterns::FileName(), "Name of the log file to be used. Old files will be overwritten!");
prm.declare_entry("Storage path", ".", Patterns::DirectoryName(), "Directory where to store the data files");
prm.declare_entry("Override files", "true", Patterns::Bool(), "If set to true, all files *tu*-files in the target directory will be erased. If not, and the directory is not empty, the program throws an error");
}
prm.leave_subsection();
prm.enter_subsection("Grid generation data");
{
prm.declare_entry("Minimal grid density", "2", Patterns::Integer(), "Value for the minimal grid density");
prm.declare_entry("Maximal grid density", "8", Patterns::Integer(), "Value for the maximum grid density");
}
prm.leave_subsection();
prm.enter_subsection("Pulse data");
{
// prm.declare_entry("Maximum pulse energy", "1e6", Patterns::Double());
// prm.declare_entry("Pulse duration", "1e-12", Patterns::Double());
// prm.declare_entry("Pulse focus size", "1e-6", Patterns::Double());
prm.declare_entry("Beam radius", "1e-3", Patterns::Double());
// prm.declare_entry("Pulse wavelength", "1.2e-6", Patterns::Double());
prm.declare_entry("Pulse calculation points", "1024", Patterns::Integer());
prm.declare_entry("Pulse wavelength", "1.2e-6", Patterns::Double());
prm.declare_entry("Focal length", "8e-3", Patterns::Double());
prm.declare_entry("Focal length inside the material", "4e-3", Patterns::Double());
prm.declare_entry("Pulse energy", "2.5e-6", Patterns::Double());
prm.declare_entry("Pulse duration", "250e-15", Patterns::Double());
}
prm.leave_subsection();
prm.enter_subsection("Simulation parameters");
{
prm.declare_entry("Absolute time step", "-1", Patterns::Double(), "If empty, a relative time step can be set");
prm.declare_entry("Relative time step", "-1", Patterns::Double(), "Time step relative to pulse length");
prm.declare_entry("Maximum inner iterations", "10", Patterns::Integer());
prm.declare_entry("Max equal values for residual", "10", Patterns::Integer());
prm.declare_entry("Theta", "0.5", Patterns::Double());
prm.declare_entry("Absolute starting time", "-1", Patterns::Double());
prm.declare_entry("Absolute maximum time", "-1", Patterns::Double());
prm.declare_entry("Relative starting time", "-1", Patterns::Double());
prm.declare_entry("Relative maximum time", "-1", Patterns::Double());
}
prm.leave_subsection();
prm.enter_subsection("Boundary parameters");
{
prm.declare_entry("Upper boundary", "1", Patterns::Integer(), "Declares boundary condition on upper boundary");
prm.declare_entry("Lower boundary", "1", Patterns::Integer(), "Declares boundary condition on lower boundary");
prm.declare_entry("Left boundary", "1", Patterns::Integer(), "Declares boundary condition on left boundary");
prm.declare_entry("Right boundary", "1", Patterns::Integer(), "Declares boundary condition on right boundary");
}
prm.leave_subsection();
}
template <int dim>
void MinimalSurfaceProblem<dim>::parse_parameters()
{
DeclExceptionMsg(ExcCreateDir, "Error while creating data directory\n");
DeclExceptionMsg(ExcDirNotEmpty, "Target directory is not empty\n");
DeclExceptionMsg(ExcOldVersion, "Version of the parameter file is too low, please update parameter file\n");
if(prm.get_integer("Version") < PARAMETER_VERSION)
AssertThrow(false, ExcOldVersion());
prm.enter_subsection("File processing data");
{
const bool use_stdcout = prm.get_bool("Use stdcout");
log_file_name = prm.get("Log file name");
(void) use_stdcout;
data_path = prm.get("Storage path");
const bool override_files = prm.get_bool("Override files");
if(data_path == ".")
data_path = "";
else
{
if((Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
{
boost::filesystem::path p(data_path);
if(!(boost::filesystem::exists(p)))
{
AssertThrow(boost::filesystem::create_directory(p) == true, ExcCreateDir());
}
else
{
if(!(boost::filesystem::is_directory(p)))
{
AssertThrow(boost::filesystem::create_directory(p) == true, ExcCreateDir());
}
else//check if directory is empty
{
if(!boost::filesystem::is_empty(p))
{
if(override_files)//Deleting all files in directory
{
for (boost::filesystem::directory_iterator end_dir_it, it(p); it!=end_dir_it; ++it) {
boost::filesystem::remove_all(it->path());
}
}
else
{
AssertThrow(false, ExcDirNotEmpty());
print_status_update(pcout, std::string("Directory is not empty, exiting\n"), true);
// throw("Directory is not empty, exiting\n");
}
}
}
}
boost::filesystem::copy_file(parameters_file_name, data_path + parameters_file_name, boost::filesystem::copy_option::overwrite_if_exists);
}
}
}
prm.leave_subsection();
prm.enter_subsection("Grid generation data");
{
min_grid_level = prm.get_integer("Minimal grid density");
max_grid_level = prm.get_integer("Maximal grid density");
}
prm.leave_subsection();
print_status_update(pcout, std::string("Got minimal and maximal grid level\n"), true);
prm.enter_subsection("Pulse data");
{
pulse_calculation_points = prm.get_integer("Pulse calculation points");
pulse_duration = prm.get_double("Pulse duration");
pulse_wavelength = prm.get_double("Pulse wavelength");
focal_length = prm.get_double("Focal length");
focal_length_inside = prm.get_double("Focal length inside the material");
pulse_energy = prm.get_double("Pulse energy");
beam_width = prm.get_double("Beam radius");
}
prm.leave_subsection();
print_status_update(pcout, std::string("Got pulse data\n"), false);
prm.enter_subsection("Simulation parameters");
{
const double abs_time_step = prm.get_double("Absolute time step");
const double rel_time_step = prm.get_double("Relative time step");
if(rel_time_step < 0 && abs_time_step < 0)
{
print_status_update(pcout, std::string("Time steps below 0, i.e. no valid values. Stopping\n"), true);
throw("Time steps below 0, i.e. no valid values. Stopping\n");
return;
}
else
if(abs_time_step < 0)
time_step = rel_time_step * pulse_duration;
else
time_step = abs_time_step;
print_status_update(pcout, std::string("Got time step values\n"), false);
max_inner_iterations = prm.get_integer("Maximum inner iterations");
max_same_res_values = prm.get_integer("Max equal values for residual");
theta = prm.get_double("Theta");
print_status_update(pcout, std::string("Got theta\n"), false);
const double abs_start_time = prm.get_double("Absolute starting time");
const double abs_max_time = prm.get_double("Absolute maximum time");
const double rel_start_time = prm.get_double("Relative starting time");
const double rel_max_time = prm.get_double("Relative maximum time");
if(abs_max_time < 0)
if(rel_max_time < 0)
{
print_status_update(pcout, std::string("Time data below 0, i.e. now valid values. Stopping\n"), true);
throw("Time data below 0, i.e. now valid values. Stopping\n");
}
else
{
start_time = rel_start_time * pulse_duration;
max_time = rel_max_time * pulse_duration;
}
else
{
start_time = abs_start_time;
max_time = abs_max_time;
}
}
prm.leave_subsection();
prm.enter_subsection("Boundary parameters");
{
boundary_conditions.push_back(prm.get_integer("Upper boundary") == 1?true:false);
boundary_conditions.push_back(prm.get_integer("Lower boundary") == 1?true:false);
boundary_conditions.push_back(prm.get_integer("Left boundary") == 1?true:false);
boundary_conditions.push_back(prm.get_integer("Right boundary") == 1?true:false);
}
prm.leave_subsection();
print_status_update(pcout, std::string("Got simulation data\n"), true);
std::stringstream time_step_stream;
time_step_stream << "Timestep is " << time_step << '\n';
print_status_update(pcout, time_step_stream.str(), true);
}
template <int dim>
void MinimalSurfaceProblem<dim>::setup_matrices(const IndexSet &partitioner, const IndexSet &relevant_partitioner)
{
system_matrix.clear();
//TrilinosWrappers::SparsityPattern sp(partitioner, partitioner, relevant_partitioner, MPI_COMM_WORLD);
//DoFTools::make_sparsity_pattern(dof_handler, sp, hanging_node_constraints, false, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD));
//DoFTools
DynamicSparsityPattern dsp(relevant_partitioner);
DoFTools::make_sparsity_pattern(dof_handler, dsp, boundary_constraints, false);
SparsityTools::distribute_sparsity_pattern(dsp, dof_handler.n_locally_owned_dofs_per_processor(),
mpi_communicator, relevant_partitioner);
system_matrix.reinit(partitioner, partitioner, dsp, mpi_communicator);
// sp.compress();
// system_matrix.reinit(sp);
}
template <int dim>
void MinimalSurfaceProblem<dim>::setup_system (const bool initial_step)
{
dof_handler.distribute_dofs (fe);
const size_t dof_numbers = dof_handler.n_dofs();
std::locale s = pcout.get_stream().getloc();
pcout.get_stream().imbue(std::locale(""));
//pcout << std::setprecision(10);
// std::string output_stream;
// output_stream = std::string("Number of active cells: ") + std::to_string(triangulation.n_global_active_cells()) + std::string(" (on ") + std::to_string(triangulation.n_levels()) + std::string(" levels)")
// + std::string("\n") + std::string("Number of degrees of freedom: ") + std::to_string(dof_numbers) + std::string("\n\n");
// print_status_update(pcout, output_stream, true);
pcout.get_stream().imbue(s);
IndexSet solution_partitioning(dof_numbers), solution_relevant_partitioning(dof_numbers);
solution_partitioning = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
//General constraint matrix
hanging_node_constraints.clear();
hanging_node_constraints.reinit(solution_relevant_partitioning);
DoFTools::make_hanging_node_constraints(dof_handler, hanging_node_constraints);
hanging_node_constraints.close();
//Newton constraint matrix;
//print_status_update(pcout, std::string("Creating newton_constraints\n"), true);
newton_constraints.clear();
newton_constraints.reinit(solution_relevant_partitioning);
DoFTools::make_hanging_node_constraints(dof_handler, newton_constraints);
VectorTools::interpolate_boundary_values(dof_handler, 0, ZeroFunction<dim>(n_components * dim), newton_constraints);
newton_constraints.close();
//Boundary constraints
//print_status_update(pcout, std::string("Creating boundary_constraints\n"), true);
boundary_constraints.clear();
boundary_constraints.reinit(solution_relevant_partitioning);
DoFTools::make_hanging_node_constraints(dof_handler, boundary_constraints);
//VectorTools::interpolate_boundary_values(dof_handler, 0, BoundaryValuesTemperatures<dim>(n_components), boundary_constraints);
VectorTools::interpolate_boundary_values(dof_handler, 0, BoundaryValuesTemperatures<dim>(n_components), boundary_constraints);
boundary_constraints.close();
setup_matrices(solution_partitioning, solution_relevant_partitioning);
system_rhs.reinit(solution_partitioning, mpi_communicator, true);
newton_update.reinit(solution_partitioning, solution_relevant_partitioning, mpi_communicator);
extension_vector.reinit(solution_partitioning, solution_relevant_partitioning, mpi_communicator);
residual.reinit(solution_partitioning, solution_relevant_partitioning, mpi_communicator);
evaluation_point.reinit(solution_partitioning, solution_relevant_partitioning, mpi_communicator);
if(initial_step)
{
present_solution.reinit(solution_partitioning, solution_relevant_partitioning, mpi_communicator);
old_solution.reinit(solution_partitioning, solution_relevant_partitioning, mpi_communicator);
}
}
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_system (void)
{
const QGauss<dim> quadrature_formula(fe.degree+1);
const FEValuesExtractors::Vector surface_TE(Electron_Temperature);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
update_gradients |
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();
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> residual_derivatives (dofs_per_cell);
for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
{
if(cell->is_locally_owned())
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
cell->get_dof_indices(local_dof_indices);
Table<3, Sacado::Fad::DFad<double> > grad_TE (n_q_points, dim, dim), grad_TE_old(n_q_points, dim, dim);
Table<2, Sacado::Fad::DFad<double> > val_TE (n_q_points, dim), val_TE_old(n_q_points, dim);
std::vector<Sacado::Fad::DFad<double> > ind_local_dof_values(dofs_per_cell);
std::vector<double> local_dof_values(cell->get_fe().dofs_per_cell);
cell->get_dof_values(present_solution, local_dof_values.begin(), local_dof_values.end());
for (size_t i = 0; i < dofs_per_cell; ++i)
{
ind_local_dof_values[i] = present_solution(local_dof_indices[i]);
ind_local_dof_values[i].diff (i, dofs_per_cell);
}
for (size_t q = 0; q < n_q_points; ++q)
for (size_t d = 0; d < dim; ++d)
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
grad_TE[q][d][local_dim] = 0;
grad_TE_old[q][d][local_dim] = 0;
}
for (size_t q = 0; q < n_q_points; ++q)
for (size_t d = 0; d < dim; ++d)
{
val_TE[q][d] = 0;
val_TE_old[q][d] = 0;
}
for (size_t q = 0; q < n_q_points; ++q)
{
for (size_t i = 0; i < dofs_per_cell; ++i)
for (size_t d = 0; d < dim; d++)
{
val_TE[q][d] += ind_local_dof_values[i] * fe_values[surface_TE].value(i, q)[d];
val_TE_old[q][d] += old_solution(local_dof_indices[i]) * fe_values[surface_TE].value(i, q)[d];
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
grad_TE[q][d][local_dim] += ind_local_dof_values[i] * fe_values[surface_TE].gradient(i, q)[d][local_dim];
grad_TE_old[q][d][local_dim] += old_solution(local_dof_indices[i]) * fe_values[surface_TE].gradient(i, q)[d][local_dim];
}
}
}
print_status_update(pcout, "Values initialized, starting assembly\n", false);
for (size_t i = 0; i < dofs_per_cell; ++i)
{
Sacado::Fad::DFad<double> R_i = 0;
for(size_t q = 0; q < n_q_points; ++q)
{
for(size_t d = 0; d < dim; ++d)
{
R_i += (val_TE[q][d] - val_TE_old[q][d])/time_step
* fe_values[surface_TE].value(i, q)[d]
* fe_values.JxW(q);
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
R_i += (theta * grad_TE[q][d][local_dim]
+ (1 - theta) * grad_TE_old[q][d][local_dim])
* fe_values[surface_TE].gradient(i, q)[d][local_dim]
* fe_values.JxW(q);
}
}
// R_i += calculate_right_side<dim, Sacado::Fad::DFad<double>>(val_TE, val_TE_old, grad_TE, grad_TE_old, fe_values[surface_TE].gradient(i, q), fe_values.JxW(q), theta);
}
for (size_t j = 0; j < dofs_per_cell; ++j)
{
residual_derivatives[j] = R_i.fastAccessDx(j);
}
//print_status_update(pcout, "Writing cell matrix\n", true);
for (size_t j = 0; j<dofs_per_cell; ++j)
cell_matrix(i, j) += residual_derivatives[j];
//print_status_update(pcout, "Setting RHS side\n", true);
cell_rhs(i) -= R_i.val();
}
cell->get_dof_indices (local_dof_indices);
newton_constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
print_status_update(pcout, std::string("Returning from assemble_matrix\n"), false);
}
// @sect4{MinimalSurfaceProblem::solve}
template <int dim>
void MinimalSurfaceProblem<dim>::solve (void)
{
IndexSet solution_relevant_partitioning(dof_handler.n_dofs());
DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
LinearAlgebraTrilinos::MPI::Vector completely_distributed_solution(dof_handler.locally_owned_dofs(), mpi_communicator);
LinearAlgebraTrilinos::MPI::Vector completely_distributed_update(dof_handler.locally_owned_dofs(), mpi_communicator);
completely_distributed_solution = present_solution;
SolverControl solver_control (dof_handler.n_dofs() * 1e1,
(system_rhs.l2_norm() > 0) ? 1e-8 * system_rhs.l2_norm() : 1e-8);
LinearAlgebraTrilinos::SolverGMRES solver (solver_control);
LinearAlgebraTrilinos::MPI::PreconditionAMG preconditioner;//Bis jetzt funktioniert IC am besten
LinearAlgebraTrilinos::MPI::PreconditionAMG::AdditionalData data;
print_status_update(pcout, "Initializing system matrix with preconditioner\n");
preconditioner.initialize(system_matrix, data);
print_status_update(pcout, "Solving\n");
solver.solve (system_matrix, completely_distributed_update, system_rhs,
preconditioner);
print_status_update(pcout, "Solving done\n");
hanging_node_constraints.distribute (completely_distributed_update);
print_status_update(pcout, "Adding to newton\n");
newton_update = completely_distributed_update;
print_status_update(pcout, std::string("L2-norm of newton-update: ") + std::to_string(completely_distributed_update.l2_norm()*1e8) + std::string("\n"), false);
newton_update.compress(VectorOperation::insert);
const double alpha = recalculate_step_length();
print_status_update(pcout, std::string("Using ") + std::to_string(alpha) + std::string(" as update\n"), true);
completely_distributed_solution.add (alpha, completely_distributed_update);
present_solution = completely_distributed_solution;
print_status_update(pcout, std::string("Returning from solve()\n"), false);
}
template <int dim>
void MinimalSurfaceProblem<dim>::refine_mesh ()
{
print_status_update(pcout, "Refining mesh\n", false);
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(fe.degree+1),
typename FunctionMap<dim>::type(),
present_solution,
estimated_error_per_cell,
ComponentMask(),
nullptr,
0,
triangulation.locally_owned_subdomain());
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.3, 0.03);
if (triangulation.n_levels() > max_grid_level)
for (auto cell = triangulation.begin_active(max_grid_level); cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();
for (auto cell = triangulation.begin_active(min_grid_level); cell != triangulation.end_active(min_grid_level); ++cell)
cell->clear_coarsen_flag ();
print_status_update(pcout, "Preparing refinement\n");
triangulation.prepare_coarsening_and_refinement ();
parallel::distributed::SolutionTransfer<dim, LinearAlgebraTrilinos::MPI::Vector> solution_transfer(dof_handler), old_solution_transfer(dof_handler);
solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
old_solution_transfer.prepare_for_coarsening_and_refinement(old_solution);
triangulation.execute_coarsening_and_refinement();
setup_system(true);
print_status_update(pcout, std::string("Interpolation\n"), false);
LinearAlgebraTrilinos::MPI::Vector distributed_solution(dof_handler.locally_owned_dofs(), mpi_communicator);
LinearAlgebraTrilinos::MPI::Vector distributed_old_solution(dof_handler.locally_owned_dofs(), mpi_communicator);
print_status_update(pcout, std::string("Doing interpolation\n"), false);
solution_transfer.interpolate(distributed_solution);
old_solution_transfer.interpolate(distributed_old_solution);
print_status_update(pcout, std::string("Adding boundary conditions\n"), false);
boundary_constraints.distribute(distributed_solution);
boundary_constraints.distribute(distributed_old_solution);
present_solution = distributed_solution;
old_solution = distributed_old_solution;
setup_system(false);
print_status_update(pcout, "Almost done with interpolation\n");
//setup_system (false);
}
// @sect4{MinimalSurfaceProblem::set_boundary_values}
// The next function ensures that the solution vector's entries respect the
// boundary values for our problem. Having refined the mesh (or just
// started computations), there might be new nodal points on the
// boundary. These have values that are simply interpolated from the
// previous mesh (or are just zero), instead of the correct boundary
// values. This is fixed up by setting all boundary nodes explicit to the
// right value:
template <int dim>
void MinimalSurfaceProblem<dim>::set_boundary_values ()
{
boundary_constraints.distribute(present_solution);
boundary_constraints.distribute(old_solution);
print_status_update(pcout, "Closed\n");
}
// @sect4{MinimalSurfaceProblem::compute_residual}
template <int dim>
double MinimalSurfaceProblem<dim>::compute_residual (const double alpha)
{
const FEValuesExtractors::Vector surface_TE(Electron_Temperature);
IndexSet solution_relevant_partitioning(dof_handler.n_dofs());
DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
evaluation_point = present_solution;
extension_vector = newton_update;
extension_vector *= alpha;
evaluation_point += extension_vector;
LinearAlgebraTrilinos::MPI::Vector local_residual(dof_handler.locally_owned_dofs(), mpi_communicator);
//evaluation_point.add (alpha, newton_update);
//evaluation_point.compress(VectorOperation::add);
print_status_update(pcout, "Creating other stuff\n");
const QGauss<dim> quadrature_formula(fe.degree+1);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
update_gradients |
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();
Vector<double> cell_residual (dofs_per_cell);
std::vector<Tensor<2, dim> > grad_TE(n_q_points), grad_TE_old(n_q_points);
std::vector<Tensor<1, dim> > val_TE(n_q_points), val_TE_old(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
print_status_update(pcout, std::string("Starting looping over cells in residual\n"), false);
for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
{
if(cell->is_locally_owned())
{
cell_residual = 0;
fe_values.reinit (cell);
fe_values[surface_TE].get_function_values(evaluation_point, val_TE);
fe_values[surface_TE].get_function_values(old_solution, val_TE_old);
fe_values[surface_TE].get_function_gradients (evaluation_point,
grad_TE);
fe_values[surface_TE].get_function_gradients (old_solution,
grad_TE_old);
print_status_update(pcout, std::string("Looping over points\n"), false);
for (size_t q = 0; q < n_q_points; ++q)
{
print_status_update(pcout, std::string("Updating residual side\n"), false);
for(size_t i = 0; i < dofs_per_cell; ++i)
{
for(size_t d = 0; d < dim; ++d)
{
cell_residual(i) -= (val_TE[q][d] - val_TE_old[q][d])/time_step
* fe_values[surface_TE].value(i, q)[d]
* fe_values.JxW(q);
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
cell_residual(i) -= (theta * grad_TE[q][d][local_dim]
+ (1 - theta) * grad_TE_old[q][d][local_dim])
* fe_values[surface_TE].gradient(i, q)[d][local_dim]
* fe_values.JxW(q);
}
}
}
}
cell->get_dof_indices (local_dof_indices);
print_status_update(pcout, std::string("Distributing values\n"), false);
boundary_constraints.distribute_local_to_global(cell_residual, local_dof_indices, local_residual);
}
}
print_status_update(pcout, std::string("Done looping\n"), false);
local_residual.compress(VectorOperation::add);
boundary_constraints.set_zero(local_residual);
local_residual.compress(VectorOperation::insert);
//std::cout << "Done looping\n";
print_status_update(pcout, std::string("Done looping II\n"), false);
//hanging_node_constraints.condense (residual);
//local_residual.compress(VectorOperation::add);
return local_residual.l2_norm();
}
template <int dim>
double MinimalSurfaceProblem<dim>::determine_step_length()
{
return alpha_val;
}
template <int dim>
double MinimalSurfaceProblem<dim>::recalculate_step_length(void)
{
std::vector<double> potential_step_sizes = {0.1, 0.25, 0.5, 0.75, 1}, potential_residual_values;
for(size_t i = 0; i < potential_step_sizes.size(); ++i)
{
auto residual = compute_residual(potential_step_sizes[i]);
potential_residual_values.push_back(residual);
}
size_t lowest_residual = std::min_element(potential_residual_values.begin(), potential_residual_values.end()) - potential_residual_values.begin();
std::stringstream lowest_residual_stream;
for(size_t i = 0; i < potential_residual_values.size(); ++i)
lowest_residual_stream << "A potential step size of " << potential_step_sizes[i] << " gives a residual of " << potential_residual_values[i] << '\n';
lowest_residual_stream << "Choosing " << potential_step_sizes[lowest_residual] << " as new alpha_val with residual_value of " << potential_residual_values[lowest_residual] << '\n';
print_status_update(pcout, lowest_residual_stream.str(), true);
return potential_step_sizes[lowest_residual];
}
template <int dim>
void MinimalSurfaceProblem<dim>::output_results (const int cycle) const
{
print_status_update(pcout, "Writing data to disk\n");
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
std::vector<std::string> solution_names = std::vector<std::string>{"solution_TE_x", "solution_TE_y"};
std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(n_components * dim, DataComponentInterpretation::component_is_scalar);
data_out.add_data_vector (old_solution, solution_names, DataOut<dim>::type_dof_data, data_component_interpretation);
Vector<float> subdomain (triangulation.n_active_cells());
for (size_t i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector (subdomain, "subdomain");
data_out.build_patches ();
const std::string filename = (data_path + "solution-" +
Utilities::int_to_string (cycle, 2) +
"." +
Utilities::int_to_string
(triangulation.locally_owned_subdomain(), 4));
std::ofstream output ((filename + ".vtu").c_str());
data_out.write_vtu (output);
if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{
std::vector<std::string> filenames;
for (size_t i = 0;
i < Utilities::MPI::n_mpi_processes(mpi_communicator);
++i)
filenames.push_back ("solution-" +
Utilities::int_to_string (cycle, 2) +
"." +
Utilities::int_to_string (i, 4) +
".vtu");
std::ofstream master_output ((data_path + "solution-" +
Utilities::int_to_string (cycle, 2) +
".pvtu").c_str());
data_out.write_pvtu_record (master_output, filenames);
}
print_status_update(pcout, "Finished writing data\n", false);
}
// @sect4{MinimalSurfaceProblem::run}
// In the run function, we build the first grid and then have the top-level
// logic for the Newton iteration. The function has two variables, one that
// indicates whether this is the first time we solve for a Newton update and
// one that indicates the refinement level of the mesh:
template <int dim>
void MinimalSurfaceProblem<dim>::run ()
{
unsigned int refinement = 0;
bool first_step = true;
GridGenerator::hyper_ball (triangulation);
static const SphericalManifold<dim> boundary;
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold (0, boundary);
//triangulation.set_boundary(1);
triangulation.refine_global(min_grid_level);
double previous_res = 0;
size_t cur_step = 0;
//bool refine_this_mesh = false;
bool next_step = false;
std::vector<double> old_residuals;
const double previous_res_threshold = 2e-3;
while (first_step || ((previous_res>previous_res_threshold) && !next_step))
{
if (first_step == true)
{
print_status_update(pcout, std::string("******** Initial mesh ********\n"), true);
print_status_update(pcout, "Setting up system\n");
setup_system (true);
print_status_update(pcout, "Setting boundary values\n");
IndexSet solution_relevant_partitioning;
LinearAlgebraTrilinos::MPI::Vector local_solution;
local_solution.reinit(dof_handler.locally_owned_dofs(), mpi_communicator);
DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
print_status_update(pcout, std::string("Before distribution\n"), false);
VectorTools::project (dof_handler,
hanging_node_constraints,
QGauss<dim>(fe.degree+1),
InitialValues<dim>(),
local_solution);
boundary_constraints.distribute(local_solution);
old_solution = local_solution;
present_solution = local_solution;
output_results(0);
//std::cout << "After distribution\n";
print_status_update(pcout, std::string("After distribution\n"), false);
first_step = false;
}
else
{
++refinement;
print_status_update(pcout, std::string("******** Refined mesh ") + std::to_string(refinement) + std::string(" ********\n"), true);
output_results(refinement * 1000 + cur_step);
refine_mesh();
}
for (size_t inner_iteration = 0; inner_iteration < max_inner_iterations; ++inner_iteration)
{
print_status_update(pcout, "Assembling system\n", false);
setup_system(false);
print_status_update(pcout, "System done\n", false);
assemble_system ();
print_status_update(pcout, "Calculating previous res\n", false);
previous_res = system_rhs.l2_norm();
print_status_update(pcout, "Trying to solve\n");
solve ();
std::stringstream L2_sstream;
L2_sstream << "Current l2 norm is " << previous_res << "\n";// << ",\tResidual (0): " << compute_residual(0) << ", Residual (0.1): " << compute_residual(0.1) << " and residual (1.0): " << compute_residual(1) << "\n";
print_status_update(pcout, L2_sstream.str(), true);
}
if(refinement > 1000)
next_step = true;
if(previous_res > previous_res_threshold)
{
if(old_residuals.size() == 0)
{
print_status_update(pcout, std::string("Added initial value ") + std::to_string(previous_res) + std::string(" to vector\n"), false);
old_residuals.push_back(previous_res);
}
else
{
if(abs(previous_res - old_residuals[old_residuals.size() - 1]) < 1e-8)
{
if(old_residuals.size() < max_same_res_values)
{
print_status_update(pcout, std::string("Adding new value ") + std::to_string(previous_res) + std::string(" to old_residuals with now a size of ") + std::to_string(old_residuals.size()) + std::string("\n"), false);
old_residuals.push_back(previous_res);
}
else
{
old_residuals.clear();
next_step = true;
}
}
else
{
print_status_update(pcout, std::string("New residual is ") + std::to_string((previous_res - old_residuals[old_residuals.size() - 1])*1e9) + std::string("\n"), false);
print_status_update(pcout, std::string("Last residual is ") + std::to_string(old_residuals[old_residuals.size() - 1]*1e9) + std::string("\n"), false);
print_status_update(pcout, std::string("Current residual is ") + std::to_string((previous_res)*1e9) + std::string("\n"), false);
old_residuals.clear();
old_residuals.push_back(previous_res);
}
}
}
}
old_solution = present_solution;
old_solution.compress(VectorOperation::insert);
output_results(1);
}
}
// @sect4{The main function}
// Finally the main function. This follows the scheme of all other main
// functions:
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step15;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
char * file_name = argv[1];
if(argc < 2)
file_name = (char*)(std::string("parameters.prm").c_str());
MinimalSurfaceProblem<2> laplace_problem_2d(file_name);
laplace_problem_2d.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;
}
parameters.prm
Description: Binary data
# Set the name of the project and target:
CMAKE_MINIMUM_REQUIRED(VERSION 3.0.0)
#C-Compilers have to be hardcoded, due to updates in CMake
set(CMAKE_C_COMPILER /usr/lib64/mpi/gcc/openmpi2/bin/mpicc CACHE PATH "" FORCE)
set(CMAKE_CXX_COMPILER /usr/lib64/mpi/gcc/openmpi2/bin/mpic++ CACHE PATH ""
FORCE)
#Get deal.II-package
find_package(deal.II 9.0.0 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
set(TARGET "main")
project(${TARGET})
cmake_policy(SET CMP0015 NEW)
cmake_policy(SET CMP0003 NEW)
# Generator specific values:
if(CMAKE_GENERATOR MATCHES "Ninja")
set(_make_command "$ ninja")
else()
set(_make_command " $ make")
endif()
#Custom targets
add_custom_target(debug
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch CMAKE_BUILD_TYPE to Debug"
)
add_custom_target(release
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch CMAKE_BUILD_TYPE to Release"
)
add_custom_target(release_production
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNO_COMMENTS")
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch CMAKE_BUILD_TYPE to Release - Production"
)
# Only mention release and debug targets if it is actually possible to
# switch between them:
if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease")
set(_switch_targets
"# ${_make_command} debug - to switch the build type to
'Debug'
# ${_make_command} release - to switch the build type to
'Release'\n"
)
ENDIF()
# Define a distclean target to remove every generated file:
add_custom_target(distclean
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target runclean
COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles
COMMAND ${CMAKE_COMMAND} -E remove
CMakeCache.txt cmake_install.cmake Makefile
build.ninja rules.ninja .ninja_deps .ninja_log
COMMENT "distclean invoked"
)
# Define a strip-comments target:
FIND_PACKAGE(Perl QUIET)
if(PERL_FOUND)
add_custom_target(strip_comments
COMMAND ${PERL_EXECUTABLE} -pi -e 's\#^[ \\t]*//.*\\n\#\#g;'
${TARGET_SRC}
COMMENT "strip comments"
)
ENDIF()
#Clean all generated files
IF("${CLEAN_UP_FILES}" STREQUAL "")
set(CLEAN_UP_FILES *.log *.*tu* *.gmv *.gnuplot *.gpl *.eps *.pov *.vtk
*.ucd *.d2 *.vtu *.pvtu)
ENDIF()
add_custom_target(runclean
COMMAND ${CMAKE_COMMAND} -E remove ${CLEAN_UP_FILES}
COMMENT "runclean invoked"
)
# Print out some usage information to file:
file(WRITE ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake
"MESSAGE(
\"###
#
# Project ${TARGET} set up with
${DEAL_II_PACKAGE_NAME}-${DEAL_II_PACKAGE_VERSION} found at
# ${DEAL_II_PATH}
#
# CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}
#
# You can now run
# ${_make_command} - to compile and link
the program
${_run_targets}#
${_switch_targets}#
")
IF(NOT CMAKE_GENERATOR MATCHES "Ninja")
FILE(APPEND
${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake
"# ${_make_command} edit_cache - to change (cached)
configuration variables
# and rerun the configure and
generate phases of CMake
#
")
ENDIF()
IF(PERL_FOUND)
FILE(APPEND
${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake
"# ${_make_command} strip_comments - to strip the source
files in this
# directory off the documentation
comments
")
ENDIF()
FILE(APPEND
${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake
"# ${_make_command} clean - to remove the
generated executable as well as
# all intermediate compilation
files
# ${_make_command} runclean - to remove all output
generated by the program
# ${_make_command} distclean - to clean the directory
from _all_ generated
# files (includes clean, runclean
and the removal
# of the generated build system)
# ${_make_command} info - to view this message
again
#
# Have a nice day!
#
###\")"
)
add_custom_target(info
COMMAND ${CMAKE_COMMAND} -P
${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake
)
# Print this message once:
IF(NOT USAGE_PRINTED)
INCLUDE(${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake)
set(USAGE_PRINTED TRUE CACHE INTERNAL "")
ELSE()
MESSAGE(STATUS "Run ${_make_command} info to print a detailed
help message")
ENDIF()
include_directories("include"
"/opt/dealii/include"
"/opt/trilinos/include"
"/opt/intel/tbb/include"
"/opt/petsc/include"
"/opt/slepc/include"
"/opt/p4est/include")
file(GLOB_RECURSE TARGET_SRC "source/*.cpp")
file(GLOB_RECURSE TARGET_INC "include/*.h" )
#get the other libraries
find_library(FILESYSTEM boost_filesystem)
find_library(MPI boost_mpi)
find_library(SERIAL boost_serialization)
find_library(ARMA armadillo)
set(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
add_executable(${TARGET} ${TARGET_SRC})
#add_dependencies (${PROJECT_NAME} pulse_propagation)
DEAL_II_SETUP_TARGET(${TARGET})
target_link_libraries(${TARGET} ${FILESYSTEM} ${MPI} ${SERIAL}
${ARMA})
