554 *
const double m = 0.5;
555 *
const double c1 = 0.;
556 *
const double c2 = 0.;
557 *
const double factor =
560 *
for (
unsigned int d = 0;
d < dim; ++
d)
571 * <a name=
"step_48-SineGordonProblemclass"></a>
615 *
const double final_time;
624 * <a name=
"step_48-SineGordonProblemSineGordonProblem"></a>
661 * <a name=
"step_48-SineGordonProblemmake_grid_and_dofs"></a>
683 *
for (
const auto &cell :
triangulation.active_cell_iterators())
684 *
if (cell->is_locally_owned())
685 *
if (cell->center().
norm() < 11)
686 *
cell->set_refine_flag();
689 *
for (
const auto &cell :
triangulation.active_cell_iterators())
690 *
if (cell->is_locally_owned())
691 *
if (cell->center().
norm() < 6)
692 *
cell->set_refine_flag();
696 *
pcout <<
" Number of global active cells: "
699 *
dof_handler.distribute_dofs(fe);
701 *
pcout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
729 *
constraints.
clear();
732 *
constraints.close();
735 *
additional_data.tasks_parallel_scheme =
738 *
matrix_free_data.reinit(mapping,
744 *
matrix_free_data.initialize_dof_vector(solution);
754 * <a name=
"step_48-SineGordonProblemoutput_results"></a>
788 *
constraints.distribute(solution);
791 *
solution.update_ghost_values();
804 *
pcout <<
" Time:" << std::setw(8) << std::setprecision(3) << time
805 *
<<
", solution norm: " << std::setprecision(5) << std::setw(7)
810 *
data_out.attach_dof_handler(dof_handler);
811 *
data_out.add_data_vector(solution,
"solution");
812 *
data_out.build_patches(mapping);
814 *
data_out.write_vtu_with_pvtu_record(
817 *
solution.zero_out_ghost_values();
824 * <a name=
"step_48-SineGordonProblemrun"></a>
825 * <
h4>SineGordonProblem::run</
h4>
851 *
pcout <<
"Number of MPI ranks: "
853 *
pcout <<
"Number of threads on each rank: "
938 *
for (time +=
time_step; time <= final_time;
945 *
wtime += timer.wall_time();
960 *
pcout <<
" Average wallclock time per time step: "
964 *
<<
"s on computations." << std::endl;
973 * <a name=
"step_48-Thecodemaincodefunction"></a>
1000 *
catch (std::exception &exc)
1002 *
std::cerr << std::endl
1004 *
<<
"----------------------------------------------------"
1006 *
std::cerr <<
"Exception on processing: " << std::endl
1007 *
<< exc.what() << std::endl
1008 *
<<
"Aborting!" << std::endl
1009 *
<<
"----------------------------------------------------"
1016 *
std::cerr << std::endl
1018 *
<<
"----------------------------------------------------"
1020 *
std::cerr <<
"Unknown exception!" << std::endl
1021 *
<<
"Aborting!" << std::endl
1022 *
<<
"----------------------------------------------------"
1045<table
align=
"center" class=
"doxtable">
1098<a name=
"step_48-Parallelrunin2Dand3D"></a><
h3>Parallel
run in 2
D and 3
D</
h3>
1110 Number
of global
active cells: 15412
1130 Time: -0.624, solution
norm: 20.488
1131 Time: -0.0381, solution
norm: 16.697
1163 Number
of global
active cells: 17592
1171 Time: -0.621, solution
norm: 123.52
1238<a name=
"step_48-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
1253 Intel's math kernel library (MKL). By using the function
1254 <code>vdSin</code> in MKL, the program uses half the computing time
1255 in 2D and 40 percent less time in 3D. On the other hand, the sine
1256 computation is structurally much more complicated than the simple
1257 arithmetic operations like additions and multiplications in the rest
1258 of the local operation.
1260 <li> <b>Higher order time stepping:</b> While the implementation allows for
1261 arbitrary order in the spatial part (by adjusting the degree of the finite
1262 element), the time stepping scheme is a standard second-order leap-frog
1263 scheme. Since solutions in wave propagation problems are usually very
1264 smooth, the error is likely dominated by the time stepping part. Of course,
1265 this could be cured by using smaller time steps (at a fixed spatial
1266 resolution), but it would be more efficient to use higher order time
1267 stepping as well. While it would be straight-forward to do so for a
1268 first-order system (use some Runge–Kutta scheme of higher order,
1269 probably combined with adaptive time step selection like the <a
1270 href="http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method">Dormand–Prince
1271 method</a>), it is more challenging for the second-order formulation. At
1272 least in the finite difference community, people usually use the PDE to find
1273 spatial correction terms that improve the temporal error.
1278<a name="step_48-PlainProg"></a>
1279<h1> The plain program</h1>
1280@include "step-48.cc"
static unsigned int n_threads()
static constexpr unsigned int rank
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static constexpr std::size_t size()
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
std::vector< index_type > data
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
constexpr types::blas_int one
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string get_current_vectorization_level()
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > cosh(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation