Commit 2ba98b61 authored by Turnhout, M.C. van's avatar Turnhout, M.C. van
Browse files

update engines matrices and equations

parent c040b4dc
......@@ -64,7 +64,7 @@ elseif istat == 2
fprintf('** Start of time step %i\n', itime);
rhsn = rhs + (M/dt - (1 - theta)*(C + K))*soln;
sol(:, itime + 1) = solvestat(q, rhsn, bndcon, [], dest);
soln = sol(:, itime+1);
soln = sol(:, itime + 1);
end
end
......
......@@ -11,23 +11,24 @@ The `engine' park (chapter \ref{introduction}, phase 2) of \mlfem{} consists of
\item read and interpret the input arrays (\kwo{coord}, \kwo{top}, \kwo{mat}, \kwo{bndcon}, etc.);
\item call the element functions to calculate the element contributions to the stiffness matrix and right-hand side of each element;
\item assemble the global stiffness matrix;
\item partition the stiffness matrix in accordance with the boundary conditions;
\item and solve the problem.
\item apply time-discretization;
\item partition the global stiffness matrix in accordance with the boundary conditions;
\item (linearize) and solve the problem.
\end{itemize}
These scripts do the actually \textsl{work}, we do not call them the `engine' without reason.\\
\subsection*{Matrices}
\section{Matrices and equations for the linear engines}
The following matrices are needed for that work and may be mentioned in the upcoming sections:
The linear engines follow the theory and terminology in \cite[chapters 14--18]{Oomens2018}, and use the following matrices and columns:
\begin{enumerate}
\item[$\mat{M}$] the (global) `mass' matrix for the final linear set of equations, for time-dependent problems
\item[$\mat{C}$] the (global) `damping' or `convective' matrix for the final linear set of equations, for first derivative terms in the PDE
\item[$\mat{K}$] the (global) `stiffness' matrix for the final linear set of equations, for second derivative terms in the PDE
\item[$\col{u}$] the (global) column with the degrees of freedom
\item[$\col{f}$] the (global) right-hand-side array for the final linear set of equations
\item[$\col{f}$] the (global) right-hand-side column for the final linear set of equations
\end{enumerate}
These global matrices are assembled from the corresponding element matrices (i.e.\ $\mat{K}$ is assembled from $\mat{K}_\mathrm{e}$, etc.) and are available for linear engines after the solution \kwo{sol} has been calculated (table \ref{enginevar}, see the element documentation in appendix \ref{elmlib} for the specifics of the element matrices for each element).
These global matrices are assembled from the corresponding element matrices (i.e.\ $\mat{K}$ is assembled from $\mat{K}_\mathrm{e}$, etc.) and are available for the user after the solution \kwo{sol} has been calculated (table \ref{enginevar}, see the element documentation in appendix \ref{elmlib} for the specifics of the element matrices for each element).
\begin{table}[h]
......@@ -38,21 +39,35 @@ These global matrices are assembled from the corresponding element matrices (i.e
engine & $\mat{M}$ & $\mat{K}$ & $\mat{C}$ & $\mat{q}$ & $\col{f}$ & $\col{u}$ & derivatives & and\dots\\
\noalign{\smallskip}\hline\noalign{\bigskip}
\kwc{fem1d} & & \kwo{K} & & &\kwo{rhs} & \kwo{sol} & \texttt{sigma} &\kwo{pos}, \kwo{dest}\\
\kwc{fem1dcd} & \kwo{M} & \kwo{K} & \kwo{C} & \kwo{q} & \kwo{rhs} & \kwo{sol}, \kwo{soln} & & \kwo{pos}, \kwo{dest}\\
\kwc{femlin\_cd} & \kwo{M} & \kwo{K} & \kwo{C} & \kwo{q} & \kwo{rhs} & \kwo{sol}, \kwo{soln} && \kwo{pos}, \kwo{dest}\\
\kwc{fem1dcd} & \kwo{M} & \kwo{K} & \kwo{C} & \kwo{q} & \kwo{rhs} & \kwo{sol} & & \kwo{pos}, \kwo{dest}\\
\kwc{femlin\_cd} & \kwo{M} & \kwo{K} & \kwo{C} & \kwo{q} & \kwo{rhs} & \kwo{sol} && \kwo{pos}, \kwo{dest}\\
\kwc{femlin\_e} & & \kwo{K} & && \kwo{rhs} & \kwo{sol} & \texttt{sigma\_elm} & \kwo{pos}, \kwo{dest}\\
\kwc{femnl} & & & & & \kwo{rhs} & \kwo{sol}, \kwo{solinc} && \kwo{pos}, \kwo{dest}\\
\kwc{femnlt} & & & && \kwo{rhs} & \kwo{sol}, \kwo{soln} && \kwo{pos}, \kwo{dest}\\
\kwc{femnlt} & & & && \kwo{rhs} & \kwo{sol} && \kwo{pos}, \kwo{dest}\\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}
The matrix $\mat{q}$ is the final (global) `stiffness' matrix that is assembled from $\mat{M}, \mat{C}$, and $\mat{K}$:
\noindent The matrix $\mat{q}$ is the final (global) `stiffness' matrix that is assembled from $\mat{M}, \mat{C}$, and $\mat{K}$ in the convection-diffusion elements:
\begin{equation}
\mat{q} = \begin{cases} \mat{C} + \mat{K} & \text{for \kwo{istat}} = 1\\[1em]
\frac{\mat{M}}{\Deltaup t}+\theta\left(\mat{C}+ \mat{K}\right) & \text{for \kwo{istat}} = 2\end{cases}
\end{equation}
So the (linear) engines solve the global linear system of equations
\begin{align}
\mat{K}\cdot \col{u} & = \col{f} &\text{\kwc{fem1d}, \kwc{femlin\_e}} \label{Kuf}\\
\mat{q}\cdot \col{u} & = \col{f} & \text{\kwc{fem1dcd}, \kwc{femlin\_cd}, \kwo{istat}} =1 \label{quf}
\end{align}
for stationary problems. The instationary convection-diffusion equations lead to \cite[eq. 15.17 \& 16.39]{Oomens2018}:
\begin{equation}
\mat{M}\pderiv{\col{u}}{t}+\left(\mat{C}+\mat{K}\right)\cdot\col{u} = \col{f} \label{prethetamat}
\end{equation}
and the convection-diffusion engines calculate the solution $u_{n+1}$ at the current time point for time-dependent problems with \cite[eq. 15.28]{Oomens2018}:
\begin{align}
\left(\frac{\mat{M}}{\Deltaup t}+\theta\left(\mat{C}+ \mat{K}\right)\right)\cdot \col{u}_{n+1} & = \left(\frac{\mat{M}}{\Deltaup t} - (1-\theta)\left(\mat{C}+\mat{K}\right)\right)\cdot \col{u}_n + \col{f} \label{thetamat}\\
\mat{q}\cdot\col{u}_{n+1} & = \col{f}_n & \text{\kwo{istat}} =2 \label{qunfn}
\end{align}
\section{1D diffusion problems: \kwc{fem1d}}\label{enginefem1d}
......@@ -63,7 +78,7 @@ The script \kwc{fem1d} uses the element \kwc{elm1d} to solve the 1D diffusion eq
\end{equation}
with the final linear set of equations \cite[equation 14.61]{Oomens2018}:
\begin{equation}
\mat{K}\cdot\col{u} = \col{f}
\mat{K}\cdot\col{u} = \col{f} \tag{\ref{Kuf}}
\end{equation}
~\\ At least the parameters $c$ and $f$, and the order of the element for \kwc{elm1d} have to be provided in \kwo{mat.mat} (section \ref{elm1d}):
......@@ -82,12 +97,13 @@ The script \kwc{fem1dcd} uses the element \kwc{elm1dcd} to solve the 1D (instati
\begin{equation}
\pderiv{u}{t} + v \pderiv{u}{x} = \pderiv{}{x}\left(c\pderiv{u}{x}\right) + f \tag{\ref{elcdeq}}
\end{equation}
with the final linear set of equations \cite[equation 15.27]{Oomens2018}:
with the final linear set of equations \cite[equation 15.27]{Oomens2018}:
\begin{align}
\mat{M}\pderiv{\col{u}}{t}+\left(\mat{C}+\mat{K}\right)\cdot\col{u} & = \col{f}\\
\mat{q}\cdot\col{u} & = \col{f}\\
\mat{M}\pderiv{\col{u}}{t}+\left(\mat{C}+\mat{K}\right)\cdot\col{u} & = \col{f} \tag{\ref{prethetamat}}\\
\mat{q}\cdot\col{u} & = \col{f} &\text{\kwo{istat}} = 1 \tag{\ref{quf}}\\
\mat{q}\cdot\col{u}_{n+1} & = \col{f}_n & \text{\kwo{istat}} = 2 \tag{\ref{thetamat}}
\end{align}
Where the $\pderiv{\col{u}}{t}$ terms only apply for the instationary case, of course.\\
\noindent At least the parameters $c$ , $v$ and $f$, and the order of the element for \kwc{elm1dcd} have to be provided in \kwo{mat.mat} (section \ref{elm1dcd}):
\begin{lstlisting}[language=mlfem,numbers=none]
......@@ -107,18 +123,18 @@ This engine is used in \kwc{demo\_1d\_diffusion.m} (the demo in chapter \ref{dem
The script \kwc{femlin\_cd} uses the element \kwc{elcd} to solve the (instationary) convection-diffusion equation (in 2D):
\begin{equation}
\pderiv{u}{t} + v \cdot\grad u = \grad\cdot\left(c\grad u\right) \tag{\ref{elcdeq}}
\pderiv{u}{t} + \vec{v} \cdot\grad u = \grad\cdot\left(c\grad u\right) \tag{\ref{elcdeq}}
\end{equation}
with the final linear set of equations \cite[equation 16.39]{Oomens2018}:
\begin{align}
\mat{M}\pderiv{\col{u}}{t}+\left(\mat{C}+\mat{K}\right)\cdot\col{u} & = \col{f}\\
\mat{q}\cdot\col{u} & = \col{f}\\
\mat{M}\pderiv{\col{u}}{t}+\left(\mat{C}+\mat{K}\right)\cdot\col{u} & = \col{f} \tag{\ref{prethetamat}}\\
\mat{q}\cdot\col{u} & = \col{f} &\text{\kwo{istat}} = 1 \tag{\ref{quf}}\\
\mat{q}\cdot\col{u}_{n+1} & = \col{f}_n & \text{\kwo{istat}} = 2 \tag{\ref{thetamat}}
\end{align}
Where the $\pderiv{\col{u}}{t}$ terms only apply for the instationary case, of course.\\
\noindent The engine uses the function \kwc{sysbuild} (section \ref{sysbuildextra}) for calculations and sets \kwo{options} $=\begin{bmatrix} 1& 1& 1 &1 & 1\end{bmatrix}$ (which should not be changed).\\
\noindent At least the parameters $c$ , $\vec{v}$ and $f$, and a parameter for \kwc{elcd} related to the topology of the element have to be provided in \kwo{mat.mat} (section \ref{elcd}):
\noindent At least the parameters $c$ and $\vec{v}$, and a parameter for \kwc{elcd} related to the topology of the element have to be provided in \kwo{mat.mat} (section \ref{elcd}):
\begin{lstlisting}[language=mlfem,numbers=none]
mat.mat(1) = c % diffusion coefficient
mat.mat(2) = vx % parameter used in elcd_a to
......@@ -149,7 +165,7 @@ The script \kwc{femlin\_e} can be used for linear problems with solids (small de
This engine solves the simple final linear set of equations
\begin{equation}
\mat{K}\cdot\col{u} = \col{f}
\mat{K}\cdot\col{u} = \col{f} \tag{\ref{Kuf}}
\end{equation}
but the structure and interpretation of $\mat{K}$ and $\col{u}$ can vary greatly between different elements (again: see the element documentation for details, appendix \ref{elmlib}).\\
......
Markdown is supported
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