Skip to content
Snippets Groups Projects
Commit 53896b0e authored by René Heß's avatar René Heß
Browse files

Cleanup exercise 05 and remove alugrid

parent 330ec39e
No related branches found
No related tags found
1 merge request!30Feature/cleanup exercises
Hide whitespace changes
Inline Side-by-side
tutorial05/exercise/doc/exercise05.tex
+ 133
− 107
View file @ 53896b0e
......@@ -40,130 +40,156 @@
\exerciseheader
\begin{Exercise}{Playing with Parameters of Adaptive the FE Algorithm}
The adaptive finite element algorithm is controlled by several parameters
which can be set in the ini-file.
Besides the polynomial degree these are:
\begin{itemize}
\item \lstinline{steps} : how many times the solve-- adapt cycle is executed.
\item \lstinline{uniformlevel} : how many times uniform mesh refinement is used before
adaptive refinement starts.
\item \lstinline{tol} : tolerance where the adaptive algorithm is stopped.
\item \lstinline{fraction} : Remember that the global error is estimated from element-wise
contributions: $$\gamma^2(u_h) = \sum_{T\in \mathcal{T}_h} \gamma_T^2(u_h).$$
The fraction parameter controls which elements are marked for refinement. More precisely,
the largest threshold $\gamma^\ast$ is determined such that
$$\sum_{\{T\in\mathcal{T}_h \,:\, \gamma_T(u_h)\geq\gamma^\ast\}} \gamma_T^2(u_h)
\geq \text{\lstinline{fraction}} \cdot \gamma^2(u_h). $$
If \lstinline{fraction} is set to zero no elements need to be refined, i.e. $\gamma^\ast=\infty$.
When \lstinline{fraction} is set to one then all elements need to be refined, i.e. $\gamma^\ast=0$.
A smaller \lstinline{fraction} parameter results in a more optimized mesh but needs more steps
and thus might be more expensive. If \lstinline{fraction} is chosen too large then almost
all elements are refined which might also be expensive, so there exists an (unknown) optimal
value for \lstinline{fraction}.
\end{itemize}
Carry out the following experiments:
\begin{enumerate}
\item Choose polynomial degree 1 and fix a tolerance, e.g. \lstinline{tol}=0.01.
Set \lstinline{steps} to a large number so that the algorithm is not stopped before the
tolerance is reached. Set \lstinline{uniformlevel} to zero. Now vary the \lstinline{fraction} parameter
and measure the execution time with
\begin{lstlisting}[basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg}]
\begin{Exercise}{Playing with the Parameters}
The adaptive finite element algorithm is controlled by several
parameters which can be set in the ini-file. Besides the polynomial
degree these are:
\begin{itemize}
\item \lstinline{steps} : how many times the solve-- adapt cycle is
executed.
\item \lstinline{uniformlevel} : how many times uniform mesh
refinement is used before adaptive refinement starts.
\item \lstinline{tol} : tolerance where the adaptive algorithm is
stopped.
\item \lstinline{fraction} : Remember that the global error is
estimated from element-wise contributions:
$$\gamma^2(u_h) = \sum_{T\in \mathcal{T}_h} \gamma_T^2(u_h).$$ The
fraction parameter controls which elements are marked for
refinement. More precisely, the largest threshold $\gamma^\ast$ is
determined such that
$$\sum_{\{T\in\mathcal{T}_h \,:\, \gamma_T(u_h)\geq\gamma^\ast\}}
\gamma_T^2(u_h) \geq \text{\lstinline{fraction}} \cdot
\gamma^2(u_h). $$
If \lstinline{fraction} is set to zero no elements need to be
refined, i.e. $\gamma^\ast=\infty$. When \lstinline{fraction} is
set to one then all elements need to be refined,
i.e. $\gamma^\ast=0$. A smaller \lstinline{fraction} parameter
results in a more optimized mesh but needs more steps and thus
might be more expensive. If \lstinline{fraction} is chosen too
large then almost all elements are refined which might also be
expensive, so there exists an (unknown) optimal value for
\lstinline{fraction}.
\end{itemize}
Carry out the following experiments:
\begin{enumerate}
\item Choose polynomial degree 1 and fix a tolerance,
e.g. \lstinline{tol}=0.01. Set \lstinline{steps} to a large
number so that the algorithm is not stopped before the tolerance
is reached. Set \lstinline{uniformlevel} to zero. Now vary the
\lstinline{fraction} parameter and measure the execution time with
\begin{lstlisting}[basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg}]
$ time ./exercise05
\end{lstlisting}
\item Repeat the previous experiment with polynomial degrees two and three. Compare the optimal
meshes generated for the different polynomial degrees.
\item The idea of the heuristic refinement algorithm refining only the elements with the
largest contribution to the error is to \textit{equilibrate} the error. To check how well this
is achieved use the calculator filter in ParaView to plot $\log(\gamma_T(u_h))$. Note that
the cell data exported by the (unchanged) code is $\gamma_T^2(u_h)$! Now check
the minimum and maximum error contributions. What happens directly at the singularity,
i.e. the point $(0,0)$?
\end{enumerate}
\end{lstlisting}%$
\item Repeat the previous experiment with polynomial degrees two and
three. Compare the optimal meshes generated for the different
polynomial degrees.
\item The idea of the heuristic refinement algorithm refining only
the elements with the largest contribution to the error is to
\textit{equilibrate} the error. To check how well this is achieved
use the calculator filter in ParaView (use cell data) to plot
$\log(\gamma_T(u_h))$. Note that the cell data exported by the
(unchanged) code is $\gamma_T^2(u_h)$! Now check the minimum and
maximum error contributions. What happens directly at the
singularity, i.e. the point $(0,0)$?
\end{enumerate}
\end{Exercise}
\begin{Exercise}{Compute the True $L_2$-error}
\textit{Disclaimer: In this exercise we compute and evalute the $L_2$-norm
of the error. The residual based error estimator implemented in the
code, however, estimates the $H^1$-seminorm (which is equivalent to the $H^1$-norm here)
and thus also optimizes the mesh with respect to that norm.
This is done to make the exercise as simple as possible. It would be more appropriate
to carry out the experiments either with the true $H^1$-norm or to change the estimator
to the $L_2$-norm.}
To do this exercise we need some more background on PDELab grid functions.
The following code known from several tutorials by now
constructs a PDELab grid function object \lstinline{g}:
\begin{lstlisting}[basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg}]
\begin{Exercise}{Compute the True $L_2$-error}
\textit{Disclaimer: In this exercise we compute and evalute the
$L_2$-norm of the error. The residual based error estimator
implemented in the code, however, estimates the $H^1$-seminorm
(which is equivalent to the $H^1$-norm here) and thus also
optimizes the mesh with respect to that norm. This is done to
make the exercise as simple as possible. It would be more
appropriate to carry out the experiments either with the true
$H^1$-norm or to change the estimator to the $L_2$-norm.}
To do this exercise we need some more background on PDELab grid
functions. The following code known from several tutorials by now
constructs a PDELab grid function object \lstinline{g}:
\begin{lstlisting}[basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg}]
Problem<RF> problem(eta);
auto glambda =
[&](const auto& e, const auto& x){return problem.g(e,x);};
[&](const auto& e, const auto& x){return problem.g(e,x);};
auto g = Dune::PDELab::makeGridFunctionFromCallable(gv,glambda);
\end{lstlisting}
Also the following code segment used in many tutorials constructs
a PDELab grid function object \lstinline{zdgf} from a grid function space and
a coefficient vector:
\begin{lstlisting}[basicstyle=\ttfamily\small,
\end{lstlisting}
Also the following code segment used in many tutorials constructs
a PDELab grid function object \lstinline{zdgf} from a grid
function space and a coefficient vector:
\begin{lstlisting}[basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg}]
typedef Dune::PDELab::DiscreteGridFunction<GFS,Z> ZDGF;
ZDGF zdgf(gfs,z);
\end{lstlisting}
PDELab grid functions can be evaluated in local coordinates given by an
element and a local coordinate within the corresponding reference element.
The result is stored in an additional argument passed by reference.
Here is a code segment doing the job:
\begin{lstlisting}[basicstyle=\ttfamily\small,
\end{lstlisting}
PDELab grid functions can be evaluated in local coordinates given
by an element and a local coordinate within the corresponding
reference element. The result is stored in an additional argument
passed by reference. Here is a code segment doing the job:
\begin{lstlisting}[basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg}]
Dune::FieldVector<RF,1> truesolution(0.0);
g.evaluate(e,x,truesolution);
\end{lstlisting}
Here we assume that \lstinline{e} is an element and \lstinline{x} is a local coordinate.
The result is stored in a \lstinline{Dune::FieldVector} with one component of type
\lstinline{RF}. This allows also to return vector-valued results.
Now here are your tasks:
\begin{enumerate}
\item Extend the code in the file \lstinline{driver.hh} in \lstinline{tutorial05/exercise/src}
to provide a grid function that provides the error $u-u_h$. Note that the method \lstinline{g}
in the parameter class already returns the true solution (for $\eta=0$).
You can solve this task by writing a lambda function subtracting the values returned
by \lstinline{g.evaluate} and \lstinline{zdgf.evaluate}.
\item Compute the $L_2$-norm of the error, i.e. $\|u-u_h\|_0 = \sqrt{\int_\Omega
(u(x)-u_h(x))^2\,dx}$. This can be done by creating a grid function
with the squared error using \lstinline{Dune::PDELab::SqrGridFunctionAdapter}
from \lstinline{<dune/pdelab/function/sqr.hh>} and the function
\lstinline{Dune::PDELab::integrateGridFunction} contained in the header file
\lstinline{<dune/pdelab/common/functionutilities.hh>}.
Output the result to the console by writing in a single line a tag
like \lstinline{L2ERROR}, the number of degrees of freedom
using \lstinline{gfs.globalSize()} and the $L_2$-norm. This
lets you grep the results later for post-processing.
\item Investigate the $L_2$-error for different polynomial degrees. Run the code
using different polynomial degrees and also uniform and adaptive (use your optimal
\lstinline{fraction}) refinement. Use \lstinline{gnuplot} to produce graphs such as
those in Figure \ref{Fig:L2}. An example gnuplot file \lstinline{plot.gp} is given.
The graph shows that for uniform refinement polynomial degree 2 is better than
polynomial degree 1 only by a constant factor. For adaptive refinement a better
convergence rate can be achieved. The plot also shows that for adaptive refinement
the optimal convergence rate can be recovered. For $P_1$ finite elements and fully regular
solution in $H^2(\Omega)$ the a-priori estimate for uniform refinement yields
$\|u-u_h\|_0\leq C h^2$. Expressing $h$ in terms of the number of degrees of freedom
$N$ yields $h\sim N^{-1/d}$, so $\|u-u_h\|_0\leq C h^2 \leq C' N^{-1}$ for $d=2$.
The curve $N^{-1}$ is shown for comparison in the plot. Clearly, the result for $P_1$ is parallel
to that line.
\item Optional: Also write the true error to the VTK output file.
\end{enumerate}
\end{lstlisting}
Here we assume that \lstinline{e} is an element and \lstinline{x}
is a local coordinate. The result is stored in a
\lstinline{Dune::FieldVector} with one component of type
\lstinline{RF}. This allows also to return vector-valued results.
Now here are your tasks:
\begin{enumerate}
\item Extend the code in the file \lstinline{driver.hh} in
\lstinline{tutorial05/exercise/task} to provide a grid function
that provides the error $u-u_h$. Note that the method
\lstinline{g} in the parameter class already returns the true
solution (for $\eta=0$). You can solve this task by writing a
lambda function subtracting the values returned by
\lstinline{g.evaluate} and \lstinline{zdgf.evaluate}.
\item Compute the $L_2$-norm of the error, i.e.
$\|u-u_h\|_0 = \sqrt{\int_\Omega (u(x)-u_h(x))^2\,dx}$. This can
be done by creating a grid function with the squared error using
\lstinline{Dune::PDELab::SqrGridFunctionAdapter} from
\lstinline{<dune/pdelab/function/sqr.hh>} and the function
\lstinline{Dune::PDELab::integrateGridFunction} contained in the
header file
\lstinline{<dune/pdelab/common/functionutilities.hh>}. Output
the result to the console by writing in a single line a tag like
\lstinline{L2ERROR}, the number of degrees of freedom using
\lstinline{gfs.globalSize()} and the $L_2$-norm. This lets you
grep the results later for post-processing.
\item Investigate the $L_2$-error for different polynomial
degrees. Run the code using different polynomial degrees and
also uniform and adaptive (use your optimal
\lstinline{fraction}) refinement. Use \lstinline{gnuplot} to
produce graphs such as those in Figure \ref{Fig:L2}. An example
gnuplot file \lstinline{plot.gp} is given.
The graph shows that for uniform refinement polynomial degree 2
is better than polynomial degree 1 only by a constant
factor. For adaptive refinement a better convergence rate can be
achieved. The plot also shows that for adaptive refinement the
optimal convergence rate can be recovered. For $P_1$ finite
elements and fully regular solution in $H^2(\Omega)$ the
a-priori estimate for uniform refinement yields
$\|u-u_h\|_0\leq C h^2$. Expressing $h$ in terms of the number
of degrees of freedom $N$ yields $h\sim N^{-1/d}$, so
$\|u-u_h\|_0\leq C h^2 \leq C' N^{-1}$ for $d=2$. The curve
$N^{-1}$ is shown for comparison in the plot. Clearly, the
result for $P_1$ is parallel to that line.
\item Optional: Also write the true error to the VTK output file.
\end{enumerate}
\end{Exercise}
\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{l2error}
......
tutorial05/exercise/solution/solution05.cc
+ 4
− 24
View file @ 53896b0e
......@@ -88,8 +88,10 @@ int main(int argc, char** argv)
ptreeparser.readOptions(argc,argv,ptree);
// read ini file
const int dim = ptree.get<int>("grid.dim");
std::string gridmanager = ptree.get("grid.manager","ug");
const int dim = 2;
std::string gridmanager = "ug";
// const int dim = ptree.get<int>("grid.dim");
// std::string gridmanager = ptree.get("grid.manager","ug");
const int degree = ptree.get<int>("fem.degree");
// use UG grid if available and selected
......@@ -113,28 +115,6 @@ int main(int argc, char** argv)
#endif
}
if (dim==2 && gridmanager=="alu")
{
#if HAVE_DUNE_ALUGRID
typedef Dune::ALUGrid<2,2,Dune::simplex,
Dune::conforming> Grid;
std::string filename = ptree.get("grid.twod.filename",
"ldomain.msh");
Dune::GridFactory<Grid> factory;
Dune::GmshReader<Grid>::read(factory,filename,true,true);
std::shared_ptr<Grid> gridp(factory.createGrid());
if (degree==1) driver<Grid,1>(*gridp,ptree);
if (degree==2) driver<Grid,2>(*gridp,ptree);
if (degree==3) driver<Grid,3>(*gridp,ptree);
if (degree==4) driver<Grid,4>(*gridp,ptree);
#else
std::cout << "You selected alugrid as grid manager "
<< "but alugrid was not found during installation"
<< std::endl;
#endif
}
if (!(gridmanager=="ug") && !(gridmanager=="alu"))
std::cout << "Example requires an unstructured grid!" << std::endl;
......
tutorial05/exercise/solution/tutorial05.ini
+ 0
− 6
View file @ 53896b0e
[grid]
dim=2
manager=ug
[grid.twod]
filename=ldomain.msh
......@@ -18,5 +14,3 @@ eta=0.0
[output]
filename=adaptive_ug
subsampling=0
tutorial05/exercise/task/exercise05.cc
+ 4
− 23
View file @ 53896b0e
......@@ -88,8 +88,10 @@ int main(int argc, char** argv)
ptreeparser.readOptions(argc,argv,ptree);
// read ini file
const int dim = ptree.get<int>("grid.dim");
std::string gridmanager = ptree.get("grid.manager","ug");
const int dim = 2;
std::string gridmanager = "ug";
// const int dim = ptree.get<int>("grid.dim");
// std::string gridmanager = ptree.get("grid.manager","ug");
const int degree = ptree.get<int>("fem.degree");
// use UG grid if available and selected
......@@ -114,27 +116,6 @@ int main(int argc, char** argv)
}
if (dim==2 && gridmanager=="alu")
{
#if HAVE_DUNE_ALUGRID
typedef Dune::ALUGrid<2,2,Dune::simplex,
Dune::conforming> Grid;
std::string filename = ptree.get("grid.twod.filename",
"ldomain.msh");
Dune::GridFactory<Grid> factory;
Dune::GmshReader<Grid>::read(factory,filename,true,true);
std::shared_ptr<Grid> gridp(factory.createGrid());
if (degree==1) driver<Grid,1>(*gridp,ptree);
if (degree==2) driver<Grid,2>(*gridp,ptree);
if (degree==3) driver<Grid,3>(*gridp,ptree);
if (degree==4) driver<Grid,4>(*gridp,ptree);
#else
std::cout << "You selected alugrid as grid manager "
<< "but alugrid was not found during installation"
<< std::endl;
#endif
}
if (!(gridmanager=="ug") && !(gridmanager=="alu"))
std::cout << "Example requires an unstructured grid!" << std::endl;
......
tutorial05/exercise/task/tutorial05.ini
+ 0
− 6
View file @ 53896b0e
[grid]
dim=2
manager=ug
[grid.twod]
filename=ldomain.msh
......@@ -18,5 +14,3 @@ eta=0.0
[output]
filename=adaptive_ug
subsampling=0
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment