Skip to content
Snippets Groups Projects
Commit a7908f2e authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on heat demo and updated the documentation.

parent ef5ea94a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 100 additions and 593 deletions
demo/init/heat.dat.1d
+ 1
− 1
View file @ a7908f2e
......@@ -42,7 +42,7 @@ heat->space->marker->ESThetaC: 0.05
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
heat->space->output->write every i-th timestep: 1
heat->space->output->write every i-th timestep: 10
heat->space->output->append index: 1
heat->space->output->index length: 6
heat->space->output->index decimals: 3
......
demo/init/heat.dat.2d
+ 1
− 1
View file @ a7908f2e
......@@ -42,7 +42,7 @@ heat->space->marker->MSGammaC: 0.1
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
heat->space->output->write every i-th timestep: 1
heat->space->output->write every i-th timestep: 10
heat->space->output->append index: 1
heat->space->output->index length: 6
heat->space->output->index decimals: 3
......
demo/init/heat.dat.3d
+ 1
− 1
View file @ a7908f2e
......@@ -43,7 +43,7 @@ heat->adapt->coarsen bisections: 3
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
heat->space->output->write every i-th timestep: 1
heat->space->output->write every i-th timestep: 10
heat->space->output->append index: 1
heat->space->output->index length: 6
heat->space->output->index decimals: 3
......
demo/src/heat.cc
+ 3
− 14
View file @ a7908f2e
......@@ -67,9 +67,7 @@ public:
void setTime(AdaptInfo *adaptInfo)
{
ProblemInstat::setTime(adaptInfo);
rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
boundaryTime = adaptInfo->getTime();
}
void closeTimestep(AdaptInfo *adaptInfo)
......@@ -120,17 +118,11 @@ public:
}
/// Returns pointer to \ref rhsTime.
double *getRHSTimePtr()
double *getRhsTimePtr()
{
return &rhsTime;
}
/// Returns pointer to \ref theta1.
double *getBoundaryTimePtr()
{
return &boundaryTime;
}
private:
/// Used for theta scheme.
double theta;
......@@ -141,9 +133,6 @@ private:
/// time for right hand side functions.
double rhsTime;
/// time for boundary functions.
double boundaryTime;
/// Pointer to boundary function. Needed for initial problem.
AbstractFunction<double, WorldVector<double> > *exactSolution;
};
......@@ -187,7 +176,7 @@ int main(int argc, char** argv)
// ===== create boundary functions =====
G *boundaryFct = new G;
boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
boundaryFct->setTimePtr(heat.getTime());
heat.setExactSolution(boundaryFct);
heatSpace.addDirichletBC(1, boundaryFct);
......@@ -195,7 +184,7 @@ int main(int argc, char** argv)
// ===== create rhs functions =====
int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
F *rhsFct = new F(degree);
rhsFct->setTimePtr(heat.getRHSTimePtr());
rhsFct->setTimePtr(heat.getRhsTimePtr());
// ===== create operators =====
......
doc/tutorial/couple.tex
+ 2
− 2
View file @ a7908f2e
......@@ -190,11 +190,11 @@ The operators for the first problem are defined like in Section
\begin{lstlisting}{}
// ===== create operators for problem2 =====
Operator matrixOperator2(problem2.getFESpace());
Operator matrixOperator2(problem2.getFeSpace());
matrixOperator2.addZeroOrderTerm(new Simple_ZOT);
problem2.addMatrixOperator(&matrixOperator2);
Operator rhsOperator2(problem2.getFESpace());
Operator rhsOperator2(problem2.getFeSpace());
rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(),
new Identity(degree)));
problem2.addVectorOperator(&rhsOperator2);
......
doc/tutorial/ellipt.tex
+ 3
− 3
View file @ a7908f2e
......@@ -202,13 +202,13 @@ The operators now are defined as follows:
\begin{lstlisting}{}
// ===== create matrix operator =====
Operator matrixOperator(ellipt.getFESpace());
Operator matrixOperator(ellipt.getFeSpace());
matrixOperator.addSecondOrderTerm(new Laplace_SOT);
ellipt.addMatrixOperator(&matrixOperator);
// ===== create rhs operator =====
int degree = ellipt.getFESpace()->getBasisFcts()->getDegree();
Operator rhsOperator(ellipt.getFESpace());
int degree = ellipt.getFeSpace()->getBasisFcts()->getDegree();
Operator rhsOperator(ellipt.getFeSpace());
rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(&rhsOperator);
\end{lstlisting}
......
doc/tutorial/heat.tex
+ 83
− 74
View file @ a7908f2e
......@@ -38,13 +38,6 @@ to the right hand side, the equation reads
\frac{u^{new}}{\tau} - \theta\Delta u^{new} = \frac{u^{old}}{\tau} + (1-\theta)\Delta u^{old} + f(\cdot, t^{old}+\theta\tau).
\end{equation}
%\begin{figure}[h]
%\center
%\includegraphics[width=0.5\textwidth]{fig/theta.pdf}
%\caption{Time discretization with the $\theta$ scheme.}
%\label{f:theta}
%\end{figure}
\begin{table}
\center
\begin{tabular}{|cc|c|c|c|}
......@@ -64,7 +57,17 @@ to the right hand side, the equation reads
\end{table}
\subsection{Source code}
Now, we describe the crucial parts of the source code. First, the functions $f$ and $g$ are defined. In contrast to the ellipt example, the functions now are time dependent. This is implemented by deriving the function classes also from class \verb+TimedObject+. This class provides a pointer to the current time, as well as corresponding setting and getting methods. The usage of a pointer to a real value allows to manage the current time in one location. All objects that deal with the same time, point to the same value. In our example, $f$ is evaluated at $t=t^{old}+\theta\tau$, while $g$ (the Dirichlet boundary function for $u^{new}$) is evaluated at $t=t^{new}$. Function $g$ is implemented as follows:
Now, we describe the crucial parts of the source code. First, the
functions $f$ and $g$ are defined. In contrast to the ellipt example,
the functions now are time dependent. This is implemented by deriving
the function classes also from class \verb+TimedObject+. This class
provides a pointer to the current time, as well as corresponding
setting and getting methods. The usage of a pointer to a real value
allows to manage the current time in one location. All objects that
deal with the same time, point to the same value. In our example, $f$
is evaluated at $t=t^{old}+\theta\tau$, while $g$ (the Dirichlet
boundary function for $u^{new}$) is evaluated at $t=t^{new}$. Function
$g$ is implemented as follows:
\begin{lstlisting}{}
class G : public AbstractFunction<double, WorldVector<double> >,
......@@ -78,7 +81,10 @@ public:
};
\end{lstlisting}
The variable $timePtr$ is a base class member of \verb+TimedObject+. This pointer has to be set once before $g$ is evaluated the first time. Implementation of function $f$ is done in the same way.
The variable $timePtr$ is a base class member of
\verb+TimedObject+. This pointer has to be set once before $g$ is
evaluated the first time. Implementation of function $f$ is done in
the same way.
\begin{figure}
\center
......@@ -86,11 +92,26 @@ The variable $timePtr$ is a base class member of \verb+TimedObject+. This pointe
\caption{UML diagram for class Heat.}
\label{f:heat_uml}
\end{figure}
Now, we begin with the implementation of class \verb+Heat+, that represents the instationary problem. In Figure \ref{f:heat_uml}, its class diagram is shown. \verb+Heat+ is derived from class \verb+ProblemInstatScal+ which leads to following properties:
Now, we begin with the implementation of class \verb+Heat+, that
represents the instationary problem. In Figure \ref{f:heat_uml}, its
class diagram is shown. \verb+Heat+ is derived from class
\verb+ProblemInstatScal+ which leads to following properties:
\begin{description}
\item \verb+Heat+ implements the \verb+ProblemTimeInterface+, so the adaptation loop can set the current time and schedule timesteps.
\item \verb+Heat+ implements \verb+ProblemStatBase+ in the role as initial (stationary) problem. The adaptation loop can compute the initial solution through this interface. The single iteration steps can be overloaded by sub classes of \verb+ProblemInstatScal+. Actually, the initial solution is computed through the method \verb+solveInitialProblem+ of \verb+ProblemTimeInterface+. But this method is implemented by \verb+ProblemInstatScal+ interpreting itself as initial stationary problem.
\item \verb+Heat+ knows another implementation of \verb+ProblemStatBase+: This other implementation represents a stationary problem which is solved within each timestep.\end{description}
\item \verb+Heat+ implements the \verb+ProblemTimeInterface+, so the
adaptation loop can set the current time and schedule timesteps.
\item \verb+Heat+ implements \verb+ProblemStatBase+ in the role as
initial (stationary) problem. The adaptation loop can compute the
initial solution through this interface. The single iteration steps
can be overloaded by sub classes of
\verb+ProblemInstatScal+. Actually, the initial solution is computed
through the method \verb+solveInitialProblem+ of
\verb+ProblemTimeInterface+. But this method is implemented by
\verb+ProblemInstatScal+ interpreting itself as initial stationary
problem.
\item \verb+Heat+ knows another implementation of
\verb+ProblemStatBase+: This other implementation represents a
stationary problem which is solved within each
timestep.\end{description}
The first lines of class \verb+Heat+ are:
......@@ -108,15 +129,21 @@ public:
}
\end{lstlisting}
The argument \verb+heatSpace+ is a pointer to the stationary problem which is solved each timestep. It is directly handed to the base class constructor of \verb+ProblemInstatScal+. In the body of the constructor, $\theta$ is read from the parameter file and stored in a member variable. The member variable \verb+theta1+ stores the value of $\theta - 1$. A pointer to this value is used later as factor in the $\theta$-scheme.
The argument \verb+heatSpace+ is a pointer to the stationary problem
which is solved each timestep. It is directly handed to the base class
constructor of \verb+ProblemInstatScal+. In the body of the
constructor, $\theta$ is read from the parameter file and stored in a
member variable. The member variable \verb+theta1+ stores the value of
$\theta - 1$. A pointer to this value is used later as factor in the
$\theta$-scheme.
The next lines show the implementation of the time interface.
\begin{lstlisting}{}
void setTime(AdaptInfo *adaptInfo)
{
ProblemInstat::setTime(adaptInfo);
rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
boundaryTime = adaptInfo->getTime();
}
void closeTimestep(AdaptInfo *adaptInfo)
......@@ -127,10 +154,12 @@ The next lines show the implementation of the time interface.
\end{lstlisting}
The method \verb+setTime+ is called by the adaptation loop to inform
the problem about the current time. The right hand side function $f$
will be evaluated at $t^{old}+\theta\tau = t^{new} - (1-\theta)\tau$,
the Dirichlet boundary function $g$ at $t^{new}$, which is the
current time.
the problem about the current time. When this function is
reimplemented, one should always call the function in the parent
class, such that all time relavant variables in \verb+ProblemInstat+
will be updated. The right hand side function $f$ will be evaluated at
$t^{old}+\theta\tau = t^{new} - (1-\theta)\tau$, the Dirichlet
boundary function $g$ at $t^{new}$, which is the current time.
The method \verb+closeTimestep+ is called at the end of each timestep
by the adaptation loop. In the default implementation of
......@@ -190,13 +219,11 @@ Now, we define some getting functions and the private member variables:
double *getThetaPtr() { return &theta; };
double *getTheta1Ptr() { return &theta1; };
double *getRHSTimePtr() { return &rhsTime; };
double *getBoundaryTimePtr() { return &boundaryTime; };
private:
double theta;
double theta1;
double rhsTime;
double boundaryTime;
AbstractFunction<double, WorldVector<double> > *exactSolution;
};
\end{lstlisting}
......@@ -212,6 +239,7 @@ int main(int argc, char** argv)
// ===== init parameters =====
Parameters::init(false, argv[1]);
Parameters::readArgv(argc, argv);
// ===== create and init stationary problem =====
ProblemScal heatSpace("heat->space");
......@@ -257,14 +285,14 @@ The functions $f$ and $g$ are declared in the following way:
\begin{lstlisting}{}
// ===== create boundary functions =====
G *boundaryFct = new G;
boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
boundaryFct->setTimePtr(heat.getTime());
heat.setExactSolution(boundaryFct);
heatSpace.addDirichletBC(1, boundaryFct);
// ===== create rhs functions =====
int degree = heatSpace.getFESpace()->getBasisFcts()->getDegree();
int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
F *rhsFct = new F(degree);
rhsFct->setTimePtr(heat.getRHSTimePtr());
rhsFct->setTimePtr(heat.getRhsTimePtr());
\end{lstlisting}
The functions interpreted as \verb+TimedObject+s are linked with the
corresponding time pointers by \verb+setTimePtr+. The boundary
......@@ -279,7 +307,7 @@ Now, we define the operators:
double zero = 0.0;
// create laplace
Operator A(heatSpace.getFESpace());
Operator A(heatSpace.getFeSpace());
A.addSecondOrderTerm(new Laplace_SOT);
A.setUhOld(heat.getOldSolution());
if (*(heat.getThetaPtr()) != 0.0)
......@@ -301,7 +329,7 @@ operator by \verb+setUhOld+.
\begin{lstlisting}{}
// create zero order operator
Operator C(heatSpace.getFESpace());
Operator C(heatSpace.getFeSpace());
C.addZeroOrderTerm(new Simple_ZOT);
C.setUhOld(heat.getOldSolution());
heatSpace.addMatrixOperator(C, heat.getInvTau(), heat.getInvTau());
......@@ -320,7 +348,7 @@ and the adaptation loop is started:
\begin{lstlisting}{}
// create RHS operator
Operator F(heatSpace.getFESpace());
Operator F(heatSpace.getFeSpace());
F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
heatSpace.addVectorOperator(F);
......@@ -356,39 +384,20 @@ heat->adapt->start time: 0.0
heat->adapt->end time: 1.0
\end{lstlisting}
Now, tolerances are determined:
Now, tolerances for the space and the time error are determined:
\begin{lstlisting}{}
heat->adapt->tolerance: 0.01
heat->adapt->rel space error: 0.5
heat->adapt->rel time error: 0.5
heat->adapt->time theta 1: 1.0
heat->adapt->time theta 2: 0.3
heat->adapt->tolerance: 0.05
heat->adapt->time tolerance: 0.05
\end{lstlisting}
The total tolerance is divided in a space tolerance and a time
tolerance. The space tolerance is the maximal allowed space error,
given by the product of \verb+tolerance+ and \verb+rel space error+. It is reached by adaptive mesh refinements. The time tolerance
is the maximal allowed error, due to the timestep size. It is given by
the product of \verb+tolerance+ and \verb+rel time error+ and
\verb+time theta 1+. It is relevant, only if an implicit time strategy
with adaptive timestep size is used. The parameter \verb+time theta 2+
is used to enlarge the timestep, if the estimated time error falls
beneath a given threshold.
\begin{lstlisting}{}
heat->adapt->strategy: 1
heat->adapt->time delta 1: 0.7071
heat->adapt->time delta 2: 1.4142
\end{lstlisting}
If \verb+strategy+ is $0$, an explicit time strategy with fixed
timestep size is used. A value of $1$ stands for the implicit
strategy. The time tolerance is reached by successively multiplying
the timestep with \verb+time delta 1+. If the estimated timestep error
is smaller than the product of \verb+tolerance+ and
\verb+rel time error+ and \verb+time theta 2+ at the end of a
timestep, the timestep size is multiplied by \verb+time delta 2+.
strategy.
The following lines determine, whether coarsening is allowed in
regions with sufficient small errors, and how many refinements or
......@@ -402,33 +411,34 @@ heat->adapt->coarsen bisections: 2
Now, the output behavior is determined:
\begin{lstlisting}{}
heat->space->output->filename: output/heat
heat->space->output->AMDiS format: 1
heat->space->output->AMDiS mesh ext: .mesh
heat->space->output->AMDiS data ext: .dat
heat->space->output->write every i-th timestep: 10
heat->space->output->append index: 1
heat->space->output->index length: 6
heat->space->output->index decimals: 3
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
heat->space->output->write every i-th timestep: 1
heat->space->output->append index: 1
heat->space->output->index length: 6
heat->space->output->index decimals: 3
\end{lstlisting}
In this example, all output filenames start with prefix
\verb+output/heat+ and end with the extensions \verb+.mesh+ and
\verb+.dat+. Output is written after every $10$th timestep. The time
of the single solution is added after the filename prefix with 6
letters, three of them are decimals. The solution for $t=0$ e.g. would
be written in the files \verb+output/heat00.000.mesh+ and
\verb+output/heat00.000.dat+.
Finally, we set parameter \verb+WAIT+ to $1$. So each call of the
macro \verb+WAIT+ in the application will lead to an interruption of
the program, until the \verb+return+ key is pressed.
\verb+output/heat+ and end with the extension \verb+.vtu+. Output is
written after every $10$th timestep. The time of the single solution
is added after the filename prefix with 6 letters, three of them are
decimals. The solution for $t=0$ e.g.\ would be written to the file
\verb+output/heat00.000.vtu+. If the parameter
\verb+ParaView animation+ is enabled, AMDiS writes for the whole
simulation one ParaView \verb+pvd+ file (in this case
\verb+output/heat.pvd+) including the names of all \verb+vtu+ files
that were created. This file makes it very easy to view and analyze
all the results of an instationary problem in ParaView.
Finally, we set parameter \verb+WAIT+ to $0$. If the variable is set
to $1$, each call of the macro \verb+WAIT+ in the application will
lead to an interruption of the program, until the \verb+return+ key is
pressed.
\begin{lstlisting}{}
WAIT: 1
WAIT: 0
\end{lstlisting}
\subsection{Macro file}
......@@ -437,8 +447,7 @@ described in Section \ref{s:ellipt macro}.
\subsection{Output}
As mentioned above, the output files look like
\verb+output/heat00.000.mesh+ and
\verb+output/heat00.000.dat+. Depending on the corresponding value in
\verb+output/heat00.000.vtu+. Depending on the corresponding value in
the parameter file only the solution after every $i$-th timestep is
written. In Figure \ref{f:heat}, the solution at three timesteps is
visualized.
......
doc/tutorial/nonlin.tex deleted 100644 → 0
+ 0
− 491
View file @ ef5ea94a
\section{Nonlinear problem}
\label{ss:nonlinear problem}
We define the nonlinear problem
\begin{eqnarray}
-\Delta u + u^4 &=& f~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
u &=& g~~~ \mbox{on } \partial\Omega.
\end{eqnarray}
We choose the functions $f$ and $g$ so that the exact solution again is $u(x)=e^{-10x^2}$. This leads to
\begin{eqnarray}
f(x)&=& -\left(400-20dow\right)e^{-10x^2} + \left(e^{-10x^2}\right)^4\\
g(x)&=&e^{-10x^2},
\end{eqnarray}
with $dow$ the world dimension.
We linearize the problem using the Newton method. First, we define an
initial guess $u_0$ of the solution which is $0$ for the first
adaptation loop iteration. In later iterations we can use the solution
of the last iteration interpolated to the current mesh as initial
guess. In each Newton step, a correction $d$ for the solution of the
last step is computed by solving
\begin{eqnarray}
\label{eq:newton step}
DF(u_n)(d) = F(u_n)
\end{eqnarray}
for $d$, where $F(u):=-\Delta u + u^4 - f$ and
\begin{eqnarray}
DF(u_n)(d) &=& \lim_{h \rightarrow 0} \frac{F(u_n+hd)-F(u_n)}{h}\\
&=& \lim_{h \rightarrow 0} \frac{-\Delta u_n - h \Delta d + \Delta u_n}{h} + \lim_{h \rightarrow 0}\frac{(u_n+hd)^4-u_n^4}{h}\\
&=& -\Delta d + 4 u_n^3 d
\end{eqnarray}
the directional derivative of $F$ at $u_n$ along $d$.
Then the solution is updated:
\begin{equation}
u_{n+1} := u_n - d.
\end{equation}
We repeat this procedure until $||d||_{L^2} < tol$ with $tol$ a given tolerance for the Newton method.
In our example, equation (\ref{eq:newton step}) reads:
\begin{equation}
-\Delta d + 4 u_n^3 d = - \Delta u_n + u_n^4 - f.
\end{equation}
\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/nonlin_iteration.pdf}
\caption{Solve step of the nonlinear problem.}
\label{f:nonlin iteration}
\end{figure}
In Figure \ref{f:nonlin iteration}, the Newton method is illustrated.
\begin{table}
\center
\begin{tabular}{|cc|c|c|c|}
\hline
& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
\hline
{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt nonlin.cc} \\
\hline
{\bf parameter file} & \tt init/ & \tt nonlin.dat.1d & \tt nonlin.dat.2d & \tt nonlin.dat.3d \\
\hline
{\bf macro file} & \tt macro/ & \tt macro\_big.stand.1d & \tt macro\_big.stand.2d & \tt macro\_big.stand.3d \\
\hline
{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt nonlin.mesh, nonlin.dat} \\
\hline
\end{tabular}
\caption{Files of the {\tt nonlin} example.}
\end{table}
\subsection{Source code}
The main idea is to realize the Newton method as implementation of the {\it ProblemIterationInterface}, which replaces the standard iteration in the adaptation loop. The Newton method has to know the problem which has to be solved, as well as the implementation of one Newton step. Estimation and adaptation is done by the problem object. The assemble step and solve step are delegated to a Newton-step object.
Now, we describe the code step by step. The function $g$ is defined like in the previous examples. In the following, we define a zero-function which is later used to implement the Dirichlet boundary condition for the Newton-step implementation (at domain boundaries no correction has to be done). The function \verb+F+ implements the right hand side function $f$.
\begin{lstlisting}{}
class Zero : public AbstractFunction<double, WorldVector<double> >
{
public:
double operator()(const WorldVector<double>& x) const
{
return 0.0;
}
};
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
/// Constructor
F(int degree)
: AbstractFunction<double, WorldVector<double> >(degree)
{}
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
int dow = x.getSize();
double r2 = x * x;
double ux = exp(-10.0 * r2);
double ux4 = ux * ux * ux * ux;
return ux4 -(400.0 * r2 - 20.0 * dow) * ux;
}
};
\end{lstlisting}
The class \verb+X3+ implements the function $u^3(x)$ used within the Newton step.
\begin{lstlisting}{}
class X3 : public AbstractFunction<double, double>
{
public:
X3() : AbstractFunction<double, double>(3) {}
/// Implementation of AbstractFunction::operator().
double operator()(const double& x) const
{
return x * x * x;
}
};
\end{lstlisting}
In the following, we define an interface which has to be implemented by the Newton-step object.
\begin{lstlisting}{}
class NewtonStepInterface
{
public:
virtual void initNewtonStep(AdaptInfo *adaptInfo) = 0;
virtual void exitNewtonStep(AdaptInfo *adaptInfo) = 0;
virtual void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) = 0;
virtual void solveNewtonStep(AdaptInfo *adaptInfo) = 0;
virtual DOFVector<double> *getCorrection() = 0;
};
\end{lstlisting}
The \verb+initNewtonStep+ method is called before each Newton step, the method \verb+exitNewtonStep+ after each Newton step. \verb+assembleNewtonStep+ assembles the linear system of equations needed for the next step, \verb+solveNewtonStep+ solves this system of equations. The solution is the correction $d$. The method \verb+getCorrection+ returns a pointer to the vector storing the correction.
Now, the Newton method will be implemented. Actually, the class \verb+NewtonMethod+ replaces the whole iteration in the adaptation loop, including mesh adaptation and error estimation. The Newton method, which is a loop over Newton steps, is one part of this iteration.
\begin{lstlisting}{}
class NewtonMethod : public ProblemIterationInterface
{
public:
NewtonMethod(const char *name,
ProblemScal *problem,
NewtonStepInterface *step)
: problemNonlin(problem),
newtonStep(step),
newtonTolerance(1e-8),
newtonMaxIter(100)
{
GET_PARAMETER(0, std::string(name) + "->tolerance", "%f",
&newtonTolerance);
GET_PARAMETER(0, std::string(name) + "->max iteration", "%d",
&newtonMaxIter);
solution = problemNonlin->getSolution();
correction = newtonStep->getCorrection();
}
\end{lstlisting}
In the constructor, pointers to the nonlinear problem and to the Newton-step object are stored to the class members \verb+problemNonlin+ and \verb+newtonStep+. Furthermore, the parameters \verb+newtonTolerance+ and \verb+newtonMaxIter+ are initialized, and pointers to the nonlinear solution and to the correction vector are stored.
The following methods define one iteration in the adaptation loop.
\begin{lstlisting}{}
void beginIteration(AdaptInfo *adaptInfo)
{
FUNCNAME("NewtonMethod::beginIteration()");
MSG("\n");
MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1);
MSG("=============================\n");
}
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
{
Flag flag = 0, markFlag = 0;
if (toDo.isSet(MARK)) markFlag = problemNonlin->markElements(adaptInfo);
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
flag = problemNonlin->refineMesh(adaptInfo);
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
flag |= problemNonlin->coarsenMesh(adaptInfo);
if (toDo.isSet(SOLVE)) {
newtonStep->initNewtonStep(adaptInfo);
int newtonIteration = 0;
double res = 0.0;
do {
newtonIteration++;
newtonStep->assembleNewtonStep(adaptInfo, flag);
newtonStep->solveNewtonStep(adaptInfo);
res = correction->L2Norm();
*solution -= *correction;
MSG("newton iteration %d: residual %f (tol: %f)\n",
newtonIteration, res, newtonTolerance);
} while ((res > newtonTolerance) && (newtonIteration < newtonMaxIter));
newtonStep->exitNewtonStep(adaptInfo);
}
if (toDo.isSet(ESTIMATE)) problemNonlin->estimate(adaptInfo);
return flag;
}
void endIteration(AdaptInfo *adaptInfo)
{
FUNCNAME("NewtonMethod::endIteration()");
MSG("\n");
MSG("end of iteration number: %d\n", adaptInfo->getSpaceIteration()+1);
MSG("=============================\n");
}
\end{lstlisting}
The methods \verb+beginIteration+ and \verb+endIteration+ only print some information to the standard output. In \verb+oneIteration+, the iteration, including the loop over the Newton steps, is defined.
Finally, the methods \verb+getNumProblems+ and \verb+getProblem+ are implemented to complete the \\\verb+ProblemIterationInterface+, and the private class members are defined.
\begin{lstlisting}{}
int getNumProblems() { return 1; };
ProblemStatBase *getProblem(int number = 0)
{
FUNCNAME("NewtonMethod::getProblem()");
if (number == 0) return problemNonlin;
ERROR_EXIT("invalid problem number\n");
return NULL;
}
private:
ProblemScal *problemNonlin;
NewtonStepInterface *newtonStep;
double newtonTolerance;
int newtonMaxIter;
DOFVector<double> *solution;
DOFVector<double> *correction;
};
\end{lstlisting}
The class \verb+Nonlin+ implements both, the nonlinear problem and the Newton-step. Since the Newton step is accessed over an own interface, it is always clear, in which role a \verb+Nonlin+ instance is called by the Newton method.
\begin{lstlisting}{}
class Nonlin : public ProblemScal,
public NewtonStepInterface
{
public:
Nonlin(const char *name)
: ProblemScal(name)
{}
\end{lstlisting}
In the constructor, the base class constructor of \verb+ProblemScal+ is called and the name is given to it.
In the initialization, the base class intialization is called, the correction vector is created and initialized, and Dirichlet boundary conditions are created.
\begin{lstlisting}{}
void initialize(Flag initFlag,
ProblemScal *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING)
{
ProblemScal::initialize(initFlag, adoptProblem, adoptFlag);
correction = new DOFVector<double>(this->getFESpace(), "old solution");
correction->set(0.0);
dirichletZero = new DirichletBC(1, &zero, feSpace_);
dirichletG = new DirichletBC(1, &g, feSpace_);
solution->getBoundaryManager()->addBoundaryCondition(dirichletG);
systemMatrix->getBoundaryManager()->
addBoundaryCondition(dirichletZero);
rhs->getBoundaryManager()->addBoundaryCondition(dirichletZero);
correction->getBoundaryManager()->addBoundaryCondition(dirichletZero);
};
\end{lstlisting}
To the solution of the nonlinear problem the function $g$ is applied as Dirichlet boundary function. The system matrix, the correction vector and the right hand side vector build the system for the Newton step. Since no correction has to be done at the domain boundaries, zero Dirichlet conditions are applied to them.
In the destructor, the allocated memory is freed.
\begin{lstlisting}{}
~Nonlin()
{
delete correction;
delete dirichletZero;
delete dirichletG;
}
\end{lstlisting}
Now, we implement the Newton step functionality. First, in \verb+initNewtonStep+, we fill the solution vector with boundary values. This will not be done automatically because we let the \verb+solution_+ pointer point to \verb+correction+. The address of \verb+solution_+ is stored in \verb+tmp+. After the Newton method is finished, the \verb+solution_+ pointer is reset to its original value in \verb+exitNewtonStep+.
\begin{lstlisting}{}
void initNewtonStep(AdaptInfo *adaptInfo)
{
solution->getBoundaryManager()->initVector(solution);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh_, -1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_BOUND);
while (elInfo) {
solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
elInfo = stack.traverseNext(elInfo);
}
solution->getBoundaryManager()->exitVector(solution);
tmp = solution;
solution = correction;
}
void exitNewtonStep(AdaptInfo *adaptInfo)
{
solution = tmp;
}
\end{lstlisting}
The implementation of \verb+assembleNewtonStep+ and \verb+solveNewtonStep+ just delegates the calls to the base class implementations in \verb+ProblemScal+.
\begin{lstlisting}{}
void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag)
{
ProblemScal::buildAfterCoarsen(adaptInfo, flag);
}
void solveNewtonStep(AdaptInfo *adaptInfo)
{
ProblemScal::solve(adaptInfo);
}
\end{lstlisting}
Finally, the \verb+getCorrection+ method is implemented and private class members are defined.
\begin{lstlisting}{}
DOFVector<double> *getCorrection()
{
return correction;
}
private:
DOFVector<double> *correction;
DOFVector<double> *tmp;
Zero zero;
G g;
DirichletBC *dirichletZero;
DirichletBC *dirichletG;
};
\end{lstlisting}
Now, we start with the main program.
\begin{lstlisting}{}
int main(int argc, char* argv[])
{
FUNCNAME("main");
TEST_EXIT(argc == 2)("usage: nonlin initfile\n");
Parameters::init(false, argv[1]);
Nonlin nonlin("nonlin");
nonlin.initialize(INIT_ALL);
AdaptInfo *adaptInfo = new AdaptInfo("nonlin->adapt", 1);
NewtonMethod newtonMethod("nonlin->newton", &nonlin, &nonlin);
AdaptStationary *adapt = new AdaptStationary("nonlin->adapt",
&newtonMethod,
adaptInfo);
\end{lstlisting}
An instance of class \verb+NewtonMethod+ was created with a pointer to
a \verb+Nonlin+ object as nonlinear problem and as Newton-step
implementation. Instead of the nonlinear problem, now, the object
\verb+newtonMethod+ is handed to the adaptation loop as implementation
of \verb+ProblemIterationInterface+.
We have to add operators representing the Newton step equation
$-\Delta d + 4u_n^3d=-\Delta u_n + u_n^4 -f$ as well as operators
representing the nonlinear problem $-\Delta u_n + u_n^4 = f$. When the
operators are given to the problem, one can determine an assemble
factor (second argument of \verb+addMatrixOperator+ and
\verb+addVectorOperator+) as well as an estimation factor (third
argument) for each operator. So, we can manage both equations in one
problem instance. Note that in the Newton step we solve for $d$. $u_n$
there is known from the last Newton step. The term $u_n^4$ was
implemented as $u_n^3 v$, where $v$ for the Newton step equation is
equal to $d$ and in the nonlinear problem equation equal to $u_n$. So,
the corresponding operator can be used in both equations just with
different factors. In Figure \ref{f:nonlin operators}, the operator
factors for the assemble and the estimate step are shown.
\begin{figure}
\center
\includegraphics[width=0.5\textwidth]{fig/nonlin_operators.pdf}
\caption{Operator factors for the assemble step and for the estimate step.}
\label{f:nonlin operators}
\end{figure}
\begin{lstlisting}{}
// ===== create operators =====
double four = 4.0;
double one = 1.0;
double zero = 0.0;
double minusOne = -1.0;
Operator *nonlinOperator0 = new Operator(nonlin.getFESpace());
nonlinOperator0->setUhOld(nonlin.getSolution());
nonlinOperator0->addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(),
new X3));
nonlin.addMatrixOperator(nonlinOperator0, &four, &one);
nonlin.addVectorOperator(nonlinOperator0, &one, &zero);
Operator *nonlinOperator2 = new Operator(nonlin.getFESpace());
nonlinOperator2->setUhOld(nonlin.getSolution());
nonlinOperator2->addSecondOrderTerm(new Laplace_SOT);
nonlin.addMatrixOperator(nonlinOperator2, &one, &one);
nonlin.addVectorOperator(nonlinOperator2, &one, &zero);
int degree = nonlin.getFESpace()->getBasisFcts()->getDegree();
Operator* rhsFunctionOperator = new Operator(nonlin.getFESpace());
rhsFunctionOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
nonlin.addVectorOperator(rhsFunctionOperator, &minusOne, &one);
\end{lstlisting}
Finally, the adaptation loop is started and after it is fnished the result is written.
\begin{lstlisting}{}
adapt->adapt();
nonlin.writeFiles(adaptInfo, true);
}
\end{lstlisting}
\subsection{Parameter file}
The used parameter file \verb+nonlin.dat.2d+ looks like:
\begin{lstlisting}{}
dimension of world: 2
nonlinMesh->macro file name: ./macro/macro_big.stand.2d
nonlinMesh->global refinements: 0
nonlin->mesh: nonlinMesh
nonlin->dim: 2
nonlin->polynomial degree: 1
nonlin->newton->tolerance: 1e-8
nonlin->newton->max iteration: 100
nonlin->solver: cg
nonlin->solver->max iteration: 1000
nonlin->solver->tolerance: 1.e-8
nonlin->solver->left precon: diag
nonlin->estimator: residual
nonlin->estimator->C0: 0.1
nonlin->estimator->C1: 0.1
nonlin->marker->strategy: 2
nonlin->marker->MSGamma: 0.5
nonlin->adapt->tolerance: 1e-1
nonlin->adapt->max iteration: 100
nonlin->output->filename: output/nonlin
nonlin->output->AMDiS format: 1
nonlin->output->AMDiS mesh ext: .mesh
nonlin->output->AMDiS data ext: .dat
\end{lstlisting}
Here, as macro file \verb+macro_big.stand.2d+ is used, which is described in Section \ref{s:nonlin macro}. The parameters \verb+nonlin->newton->tolerance+ and \verb+nonlin->newton->max iteration+ determine the tolerance of the Newton solver and the maximal number of Newton steps within each iteration of the adaptation loop.
\subsection{Macro file}
\label{s:nonlin macro}
The used macro file \verb+macro_big.stand.2d+ is very similar to the macro file of the first example described in Section \ref{s:ellipt macro}. Only the domain was changed from $\Omega=[0,1]^2$ to $\Omega=[-1, 1]^2$ by adapting the vertex coordinates correspondingly.
\subsection{Output}
\begin{figure}
\center
\begin{tabular}{cc}
\includegraphics[width=0.4\textwidth]{fig/nonlin_solution.jpg}&
\includegraphics[width=0.4\textwidth]{fig/nonlin_mesh.jpg}\\
solution & mesh
\end{tabular}
\caption{Solution and final mesh of the nonlinear problem}
\label{f:nonlin mesh and solution}
\end{figure}
In Figure \ref{f:nonlin mesh and solution}, the solution and the final mesh written after the adaptation loop are visualized. The solution is shown as height field where the values are interpreted as z-coordinates.
\ No newline at end of file
doc/tutorial/parametric.tex
+ 1
− 1
View file @ a7908f2e
......@@ -170,7 +170,7 @@ Before we let the torus rotate, some variables are defined. We set the rotation
DegreeOfFreedom dof;
WorldVector<double> x;
const FiniteElemSpace *feSpace = torus.getFESpace();
const FiniteElemSpace *feSpace = torus.getFeSpace();
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numBasFcts = basFcts->getNumber();
DegreeOfFreedom *localIndices = new double[numBasFcts];
......
doc/tutorial/tutorial.pdf
+ 0
− 0
View file @ a7908f2e
No preview for this file type
doc/tutorial/vecellipt.tex
+ 5
− 5
View file @ a7908f2e
......@@ -79,12 +79,12 @@ The operator definitions for the first equation are:
\begin{lstlisting}{}
// ===== create operators =====
Operator matrixOperator00(vecellipt.getFESpace(0), vecellipt.getFESpace(0));
Operator matrixOperator00(vecellipt.getFeSpace(0), vecellipt.getFeSpace(0));
matrixOperator00.addSecondOrderTerm(new Laplace_SOT);
vecellipt.addMatrixOperator(&matrixOperator00, 0, 0);
int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree();
Operator rhsOperator0(vecellipt.getFESpace(0));
int degree = vecellipt.getFeSpace(0)->getBasisFcts()->getDegree();
Operator rhsOperator0(vecellipt.getFeSpace(0));
rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
vecellipt.addVectorOperator(&rhsOperator0, 0);
\end{lstlisting}
......@@ -105,11 +105,11 @@ operator vector at position $0$.
Now, the operators for the second equation are defined:
\begin{lstlisting}{}
Operator matrixOperator10(vecellipt.getFESpace(1), vecellipt.getFESpace(0));
Operator matrixOperator10(vecellipt.getFeSpace(1), vecellipt.getFeSpace(0));
matrixOperator10.addZeroOrderTerm(new Simple_ZOT);
vecellipt.addMatrixOperator(matrixOperator10, 1, 0);
Operator matrixOperator11(vecellipt.getFESpace(1), vecellipt.getFESpace(1));
Operator matrixOperator11(vecellipt.getFeSpace(1), vecellipt.getFeSpace(1));
matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0));
vecellipt.addMatrixOperator(matrixOperator11, 1, 1);
\end{lstlisting}
......
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