\item Work on the problem presented in the \lstinline{tutorial04}.
\item Work on the problem presented in \lstinline{tutorial04}.
\item Try various time integrators, in particular the Crank-Nicolson method.
\item Explore polynomial degree greater than $2$ by changing the blocking
\item Explore polynomial degrees greater than $2$ by changing the blocking
to \lstinline{none}.
\item Implement your own local operator that represents elliptic projection method (see Tutorial notes and Eriksson et al. \cite{Eriksson})
\item Implement your own local operator that switches to an elliptic projection method (see the tutorial notes and Eriksson et al. \cite{Eriksson})
\item Check the energy conservation property in the numerical scheme.
\end{itemize}
...
...
@@ -58,7 +58,7 @@ to \lstinline{none}.
\lstset{language=bash}
The code of \textbf{exercise04} solves the the wave equation formulated as a first order in time system. As already explained in the \lstinline{tutorial04} by renaming $u_0=u$ and introducing $u_1=\partial_t u_0=\partial_t u$ we can write the wave equation as a system of two equations:
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$:
\begin{subequations}
\label{eq:SystemForm1}
\begin{align}
...
...
@@ -71,34 +71,34 @@ u_1 &= w &&\text{at $t=0$}.
\end{align}
\end{subequations}
Since $u_0=u=0$ on the boundary we also have $\partial_t u = u_1=0$ on the boundary.
But one may also omit the boundary condition on $u_1$.
Alternatively, one may also simply omit the boundary condition on $u_1$.
The code to this exercise can be recompiled individually \textbf{in
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/src
[user@localhost]$ make
\end{lstlisting}
The structure of the code is very similar to the previous tutorials,
it 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.cc! -- main program,
\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!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 the tutorial03. The \lstinline!pvd!-file can
be visualized by ParaView and it consists of a collection of the
corresponding \lstinline!vtu!-files. To see 1D solution, one can apply ''Plot Over Line '' filer, note that our solution is depicted by $u_0$.
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}
...
...
@@ -110,55 +110,51 @@ But one may also omit the boundary condition on $u_1$.
\begin{lstlisting}
Dune::PDELab::Alexander2Parameter<RF> pmethod;
\end{lstlisting}
Recall the previous \lstinline!exercise03! and change the one
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 "Append Attributes" filter before "Plot Over Line". Do not forget to change the torder parameter in \lstinline!exercise04.ini! in order to get the right scheme.
\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! file. 2d give you nicer pictures under the cost of longer runtime.
\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 runtimes.
\begin{Exercise}{ Explore polynomial degree greater than $2$ by changing the blocking
\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 the polynomial degree 3
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.
\textbf{Note that}\lstinline!deg! is a static template parameter.
\textit{Program should compile and give a segmentation error when started.}
\textbf{Step 2:}
Before we start please have a look at \lstinline!tutorial04! description of driver and read the part that describes specification of an ordering when creating
product spaces. We recall \textit{Important notice}: Using fixed block structure in ISTL requires that there is the same number of degrees of freedom per entity. This is true for polynomial degrees one and two but not for higher polynomial degree!
To avoid segfault 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 degree!
See how it is done for the $\Delta u_0$, we recall corresponding part of \lstinline{WaveL2}\\\lstinline{alpha_volume} method
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++) {
...
...
@@ -214,32 +210,32 @@ for (size_t i=0; i<lfsu0.size(); i++) {
}
\end{lstlisting}
\item Do not forget about proper changes in\lstinline{jacobian_volume} methods. Since the problem is linear this should not cost you too much effort.
\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 energy $$E(t)=\|u_1\|^2+\|\nabla u_0\|^2$$
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 Dune PDELab functions
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
\end{lstlisting}
If you have problems with this task check online documentation
If you have problems with this task check the online documentation at
\url{https://www.dune-project.org/doc/doxygen/html/modules.html} or the solution.
\end{Exercise}
\textbf{Additional task}\\
If you are done with this exercises, you can play with 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 proper lambda function)
\item Change the initial conditions to $\sin(2x)$ (implement it as a lambda function)
