The code of \textbf{exercise04} solves the wave equation formulated as a first order in time system. As already explained in \lstinline{tutorial04}, we can write the wave equation as a system of two equations by substituting $u_0=u$ and introducing $u_1=\partial_t u_0=\partial_t u$:
$\partial_t u = u_1=0$ on the boundary. Alternatively, one may
also simply omit the boundary condition on $u_1$.
The code forthis exercise can be recompiled individually \textbf{in
your build directory} by typing make:
\begin{lstlisting}
[user@localhost]$ cd release-build/dune-pdelab-tutorials/tutorial04/exercise/task
[user@localhost]$ make clean && make
[user@localhost]$ cd release-build/dune-pdelab-tutorials/tutorial04/exercise/task
[user@localhost]$ make clean && make
\end{lstlisting}
Note that in contrast to earlier exercises, three different executables are built from the same source files,
one each for 1D, 2D and 3D, in order to reduce compile times and compiler memory usage.
Note that in contrast to earlier exercises, three different
executables are built from the same source files, one each for 1D,
2D and 3D, in order to reduce compile times and compiler memory
usage.
The structure of the code is very similar to the previous tutorials and
consists of the following files:
The structure of the code is very similar to the previous tutorials
and consists of the following files:
\begin{itemize}
\item\lstinline!exercise04.ini! -- holds parameters read by various parts of the code which control the execution,
\item\lstinline!exercise04.ini! -- holds parameters read by various
parts of the code which control the execution,
\item\lstinline!exercise04.cc! -- main program,
\item\lstinline!driver.hh! -- instantiates the necessary PDELab classes for solving a linear instationary problem and finally solves the problem,
\item\lstinline!wavefem.hh! -- contains the local operator classes \lstinline!WaveFEM! and \lstinline!WaveL2! realizing the spatial and temporal residual forms $r(u,v,t)$ and $m(u,v,t)$, respectively.
\item\lstinline!driver.hh! -- instantiates the necessary PDELab
classes for solving a linear instationary problem and finally
solves the problem,
\item\lstinline!wavefem.hh! -- contains the local operator classes
\lstinline!WaveFEM! and \lstinline!WaveL2! realizing the spatial
and temporal residual forms $r(u,v,t)$ and $m(u,v,t)$,
respectively.
\end{itemize}
As in the previous exercises you can control most of the settings
through the ini-file \lstinline!exercise04.ini!. Get an overview of
the configurable settings, compile and run \lstinline!exercise04!.
The program writes output with the extension \lstinline!pvd!. This is one
of several ways to write VTK output for the instationary case,
c.f.~the documentation of tutorial03. The \lstinline!pvd!-file can
be visualized by ParaView and consists of a collection of the
corresponding \lstinline!vtu!-files. In order to visualize a 1D solution, one can apply the ''Plot Over Line '' filter. Note that our solution is always given by $u_0$.
The program writes output with the extension \lstinline!pvd!. This
is one of several ways to write VTK output for the instationary
case, c.f.~the documentation of tutorial03. The \lstinline!pvd!-file
can be visualized by ParaView and consists of a collection of the
corresponding \lstinline!vtu!-files. In order to visualize a 1D
solution, one can apply the ''Plot Over Line '' filter. Note that
our solution is always given by $u_0$.
\end{Exercise}
\begin{Exercise}{ Try various time integrators, in particular the Crank-Nicolson method}
\begin{Exercise}{Try various time integrators, in particular the Crank-Nicolson method}
We want to examine the numerical solution under different time
discretization schemes, eg. Implicit Euler or Crank-Nicolson. In order to change the time discretization scheme you will have to go to the file \lstinline!driver.hh! and
search for the line
discretization schemes, eg. Implicit Euler or Crank-Nicolson. In
order to change the time discretization scheme you will have to go
to the file \lstinline!driver.hh! and search for the line
\begin{lstlisting}
Dune::PDELab::Alexander2Parameter<RF> pmethod;
Dune::PDELab::Alexander2Parameter<RF> pmethod;
\end{lstlisting}
Recall the previous \lstinline!exercise03! and change the one
step $\theta$ scheme in a way that corresponds to the Crank-Nicolson method.
\textbf{Note that} you can compare the solutions in ParaView. To do that you have to rename the second solution in \lstinline!exercise04.ini! and use the "Append Attributes" filter before "Plot Over Line". Do not forget to change the parameter \lstinline!torder! in \lstinline!exercise04.ini! in order to get the correct scheme.
Recall the previous \lstinline!exercise03! and change the one step
$\theta$ scheme in a way that corresponds to the Crank-Nicolson
method.
\textbf{Note that} you can compare the solutions in ParaView. To do
that you have to rename the second solution in
\lstinline!exercise04.ini! and use the "Append Attributes" filter
before "Plot Over Line". Do not forget to change the parameter
\lstinline!torder! in \lstinline!exercise04.ini! in order to get the
correct scheme.
\end{Exercise}
\textbf{Remark} You can decide about the domain dimension, see \lstinline!exercise04.ini!. A 2D simulation will give you nicer pictures at the cost of longer run times.
\textbf{Remark} You can decide about the domain dimension, see
\lstinline!exercise04.ini!. A 2D simulation will give you nicer
pictures at the cost of longer run times.
\begin{Exercise}{Explore polynomial degrees greater than $2$ by changing the blocking
to \lstinline{none}.}
\textbf{Step 1:}
Find the place where your Local Finite Element Maps are created and make it possible to use polynomial degree 3
\begin{Exercise}{Explore polynomial degrees greater than $2$ by
changing the blocking to \lstinline{none}.}
\textbf{Step 1:} Find the place where your Local Finite Element Maps
are created and make it possible to use polynomial degree 3
\textbf{Note that}\lstinline!deg! is a static template parameter.
Compile your program and check the results.
\textit{Program should compile and give a segmentation error when started.}
\textbf{Note that}\lstinline!deg! is a static template parameter.
\textbf{Step 2:}
\textit{Program should compile and give a segmentation error when
started.}
Before we start please have a look at the description of the driver in \lstinline!tutorial04! and read the part that describes the specification of an ordering when creating
product spaces. Using fixed block structure in ISTL requires that the number of degrees of freedom per entity is constant for each geometry type. This is true for polynomial degrees one and two but not for higher polynomial degree!
\textbf{Step 2:}
To avoid segfaults you need to change
Before we start please have a look at the description of the driver
in \lstinline!tutorial04! and read the part that describes the
specification of an ordering when creating product spaces. Using
fixed block structure in ISTL requires that the number of degrees of
freedom per entity is constant for each geometry type. This is true
for polynomial degrees one and two but not for higher polynomial
See how it is done for the $\Delta u_0$, we recall the corresponding part of \lstinline{WaveL2::alpha_volume()}:
See how it is done for the $\Delta u_0$, we recall the
corresponding part of \lstinline{WaveL2::alpha_volume()}:
\begin{lstlisting}
// integrate u*phi_i
for (size_t i=0; i<lfsu0.size(); i++) {
...
...
@@ -213,30 +243,42 @@ for (size_t i=0; i<lfsu0.size(); i++) {
}
\end{lstlisting}
\item Do not forget to update the \lstinline{jacobian_volume()} methods accordingly. As the problem is linear, this should not be too difficult.
\end{itemize}
\item Do not forget to update the \lstinline{jacobian_volume()}
methods accordingly. As the problem is linear, this should not be
too difficult.
\end{itemize}
\end{Exercise}
\begin{Exercise}{Energy conservation}
If we multiply \eqref{eq:2a} by $u_1$ and \eqref{eq:elip} by $u_0$ and add, the terms $-(\Delta u_0,u_1)$ and $(\Delta u_1,u_0)$ cancel out, leading to the conclusion that the energy $$E(t)=\|u_1\|^2+\|\nabla u_0\|^2$$
is constant in time.
If we multiply \eqref{eq:2a} by $u_1$ and \eqref{eq:elip} by $u_0$
and add, the terms $-(\Delta u_0,u_1)$ and $(\Delta u_1,u_0)$
cancel out, leading to the conclusion that the energy
$$E(t)=\|u_1\|^2+\|\nabla u_0\|^2$$ is constant in time.
Your task is to check the energy conservation for the elliptic projection and Crank-Nicolson in time. You can use the following PDELab utilities:
Your task is to check the energy conservation for the elliptic
projection and Crank-Nicolson in time. You can use the following
PDELab utilities:
\begin{lstlisting}
Dune::PDELab::SqrGridFunctionAdapter
Dune::PDELab::integrateGridFunction
Dune::PDELab::DiscreteGridFunctionGradient
Dune::PDELab::SqrGridFunctionAdapter
Dune::PDELab::SqrGridFunctionAdapter
Dune::PDELab::integrateGridFunction
Dune::PDELab::DiscreteGridFunctionGradient
Dune::PDELab::SqrGridFunctionAdapter
\end{lstlisting}
If you have problems with this task check the online documentation at
\url{https://www.dune-project.org/doxygen/2.5.0/modules.html} or the solution.
If you have problems with this task check the online documentation
at \url{https://www.dune-project.org/doxygen/2.5.0/modules.html} or
the solution.
\end{Exercise}
\textbf{Additional task}\\
If you are done with these exercises, you can play with the initial conditions. Use the following setting in one space dimension:
If you are done with these exercises, you can play with the initial
conditions. Use the following setting in one space dimension:
\begin{itemize}
\item Change the initial conditions to $\sin(2x)$ (implement it as a lambda function)
\item Change the initial conditions to $\sin(2x)$ (implement it as a
