Medusa: Coordinate Free Mehless Method implementation - User contributions [en]

Solving sparse systems

2023-02-15T07:42:37Z

Mitjajancic:

There are many methods available for solving sparse systems. We compare some of them here.
<figure id="fig:matrix1">
[[File:matrix.png|300px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Matrix of the discretized PDE. </caption>]]
</figure>
Mathematica has the following methods available (https://reference.wolfram.com/language/ref/LinearSolve.html#DetailsAndOptions)
* direct: banded, cholesky, multifrontal (direct sparse LU)
* iterative: Krylov

Matlab has the following methods:
* direct: https://www.mathworks.com/help/matlab/ref/mldivide.html#bt42omx_head
* iterative: https://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html#brzoiix, including bicgstab, gmres

Eigen has the following methods: (https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
* direct: sparse LU
* iterative: bicgstab, cg

Solving a simple sparse system $A x = b$ for steady space of heat equation in 1d with $n$ nodes, results in a matrix shown in <xr id="fig:matrix1"/>.

The following timings of solvers are given in seconds:
{| class="wikitable"
|-
! $n = 10^6$
! Matlab
! Mathematica
! Eigen
|-
! Banded
| 0.16
| 0.28
| 0.04
|-
! SparseLU
| /
| 1.73
| 0.82
|-
|}

Incomplete LU preconditioner was used for BICGStab.
Without the preconditioner BICGStab does not converge.

==Parallel execution ==
BICGStab can be run in parallel, as explain in the general parallelism: https://eigen.tuxfamily.org/dox/TopicMultiThreading.html, and specifically

'''"When using sparse matrices, best performance is achied for a row-major sparse matrix format. Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled".'''

Eigen uses number of threads specified my OopenMP, unless <code>Eigen::setNbThreads(n);</code> was called.
Minimal working example:

<figure id="fig:eigen_par">
[[File:eigen_parallel_2.png|1200px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Memory and CPU usage. C stands for construction of the system, L stands for the calculation of ILUT preconditioner and S for BICGStab iteration.</caption>]]
</figure>

<syntaxhighlight lang="c++" line>
#include <iostream>
#include <vector>
#include "Eigen/Sparse"
#include "Eigen/IterativeLinearSolvers"

using namespace std;
using namespace Eigen;

int main(int argc, char* argv[]) {

assert(argc == 2 && "Second argument is size of the system.");
stringstream ss(argv[1]);
int n;
ss >> n;
cout << "n = " << n << endl;

// Fill the matrix
VectorXd b = VectorXd::Ones(n) / n / n;
b(0) = b(n-1) = 1;
SparseMatrix<double, RowMajor> A(n, n);
A.reserve(vector<int>(n, 3)); // 3 per row
for (int i = 0; i < n-1; ++i) {
A.insert(i, i) = -2;
A.insert(i, i+1) = 1;
A.insert(i+1, i) = 1;
}
A.coeffRef(0, 0) = 1;
A.coeffRef(0, 1) = 0;
A.coeffRef(n-1, n-2) = 0;
A.coeffRef(n-1, n-1) = 1;

// Solve the system
BiCGSTAB<SparseMatrix<double, RowMajor>, IncompleteLUT<double>> solver;
solver.setTolerance(1e-10);
solver.setMaxIterations(1000);
solver.compute(A);
VectorXd x = solver.solve(b);
cout << "#iterations: " << solver.iterations() << endl;
cout << "estimated error: " << solver.error() << endl;
cout << "sol: " << x.head(6).transpose() << endl;

return 0;
}
</syntaxhighlight>

was compiled with <code>g++ -o parallel_solve -O3 -fopenmp solver_test_parallel.cpp</code>.
<xr id="fig:eigen_par"/> was produced when the program above was run as <code>./parallel_solve 10000000</code>.

==ILUT factorization==
This is a method to compute a general algebraic preconditioner which has the ability to control (to some degree) the memory and time complexity of the solution.
It is described in a paper by

Saad, Yousef. "ILUT: A dual threshold incomplete LU factorization." Numerical linear algebra with applications 1.4 (1994): 387-402.

It has two parameters, tolerance $\tau$ and fill factor $f$. Elements in $L$ and $U$ factors that would be below relative tolerance $\tau$ (say, to 2-norm of this row) are dropped.
Then, only $f$ elements per row in $L$ and $U$ are allowed and the maximum $f$ (if they exist) are taken. The diagonal elements are always kept. This means that the number of non-zero
elements in the rows of the preconditioner can not exceed $2f+1$.

Greater parameter $f$ means slover but better $ILUT$ factorization. If $\tau$ is small, having a very large $f$ means getting a complete $LU$ factorization. Typical values are $5, 10, 20, 50$.
Lower $\tau$ values mean that move small elements are kept, resulting in a more precise, but longer computation. With a large $f$, the factorization approaches the direct $LU$ factorization as $\tau \to 0$. Typical values are from 1e-2
to 1e-6.

Whether the method converges or diverges for given parameters is very sudden and should be determined for each apllication separately. An example of
such a convergence graph for the diffusion equation is presented in figures below:

[[File:bicgstab_err.png|400px]][[File:bicgstab_conv.png|400px]][[File:bicgstab_iter.png|400px]]

We can see that low BiCGSTAB error estimate coincides with the range of actual small error compared to the analytical solution as well as with the range of lesser number of iterations.

==Pardiso library==
[[Pardiso]]

==Herzian contact ==
The matrix of the [[Hertzian contact]] problem is shown below:
[[File:matrix_hertzian.png|500px]]

The condition number iz arther large and grows superlinearly with size in the example below.
[[File:cond.png|600px]]

= Discussion based on practical observations =

In this section we analyse different solvers, all available as part of
the [https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html Eigen library] supported by
the [https://e6.ijs.si/medusa/ Medusa]. Namely we analyse the performance
of two direct solvers

* [https://eigen.tuxfamily.org/dox/classEigen_1_1SparseLU.html SparseLU]
* [https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html SparseQR]

three combinations of general purpose iterative
solver [https://eigen.tuxfamily.org/dox/classEigen_1_1BiCGSTAB.html BiCGSTAB]

* BiCGSTAB
* BiCGSTAB with ILUT preconditioner
* BiCGSTAB with ILUT preconditioner and exploited '''[https://eigen.tuxfamily.org/dox/classEigen_1_1SolveWithGuess.html solveWithGuess()]''' functionality

and a wrapper to external direct solver

* [https://eigen.tuxfamily.org/dox/classEigen_1_1PardisoLU.html PardisoLU].

The performance tests are focused on the accuracy of the numerical solution, number of iterations and time required to solve the final sparse system. We also comment on the amount of RAM required by a specific solver. If required, the [https://gitlab.com/e62Lab/Medusa-disucssions/-/tree/master/eigen_solvers gitlab repository] with the implementation is also provided.

<syntaxhighlight lang="text">
NOTE: Although the SparseQR has been included in our tests, significantly larger solve times compared to other solvers have been observed. For that reason, the SparseQR is discarded from further analyses even though it was able to obtain a solution of comparable accuracy.
</syntaxhighlight>

<syntaxhighlight lang="text">
Computational times: To obtain the computational times the code was run on a single CPU with assured CPU affinity using the ''taskset'' command.
</syntaxhighlight>

== The two-dimensional Poisson problem ==

First, we solve the Poisson problem
$$\nabla ^2 u(\boldsymbol x) = f(\boldsymbol x),$$

for $\boldsymbol x \in \Omega $ where the domain $\Omega$ is a unit $d=2$ dimensional circle with radius $R=1$.
We use mixed boundary conditions, i.e., Dirichlet boundary for all $x < 0$ and Neumann boundary for all $x \geq 0$. Additionally, ghost or fictitious nodes were added to the Neumann boundary.

We perform analyses on two scenarios, i.e., two definitions of $f(\boldsymbol x)$:

* $f(\boldsymbol x) = -d\pi ^2 \prod_{i = 1}^d \sin(\pi x_i)$ yielding
** $f_d(\boldsymbol x) = \prod_{i = 1}^d \sin(\pi x_i)$ on the Dirichlet boundary and
** $ \boldsymbol f_n(\boldsymbol x) = \pi \sum_{i = 1}^d \cos(\pi x_i) \prod_{j=1, j\neq i}^d \sin (\pi x_j)
\boldsymbol e_i$ with unit vectors $\boldsymbol e_i$ on the Neumann boundary.
* $f(\boldsymbol x) = 2a e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}(2a\left \| \boldsymbol x -
\boldsymbol x_s \right \| - d)$ for width $a=10^2$ of the exponential source positioned at $\boldsymbol x_s=\boldsymbol 0$ yielding
** $f_d(\boldsymbol x) = e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$ yielding on Dirichlet boundary and
** $\boldsymbol f_n (\boldsymbol x) = -2a(\boldsymbol x - \boldsymbol x_s)e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$ on Neumann boundary.

Example numerical solutions to both scenarios are shown in figure below.

[[File:solution_poisson.png|800px]]

=== Solver performance ===

Below, we show the performance of the selected solvers in case of a two-dimensional Poisson problem. The first column shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of iterations a given solver required (direct solvers were assigned a value 1) and the third column shows the total time required to solve the system and obtain the numerical solution $\widehat u$.

[[File:convergence_poisson_sin.png|1200px]]

[[File:convergence_poisson_exp.png|1200px]]

In terms of accuracy, all solvers under consideration perform equally good, with an exception of BiCGSTAB with ILUT preconditioner, who failed to obtain a solution of sufficient accuracy at the finest domain discretization with approximately $2\cdot 10^6$ nodes. The reason is most likely in the ILUT preconditioner, set to limit solver iterations to 200 with fill factor 40, drop tolerance $10^{-5}$ and global tolerance $10^{-15}$. A different ILUT setup would most probably yield a more desirable behaviour, but it nevertheless shows that when using a ILUT preconditioner one must be careful before making conclusions. It is also important to note that the SparseLU ran out of the available 11.5 Gb RAM and, therefore, failed to obtain a numerical solution. Other solvers had no issues with the available RAM space.

It is perhaps more interesting to study the number of iterations required by the solvers. Obviously, direct solvers require a "single iteration". Iterative solvers, on the other hand, require multiple iterations. Generally speaking, we observe that the number of solver iterations increases with the number of discretization points. However, it is interesting to see that the pure BiCGSTAB shows the most stable and predictable behaviour. Despite the fact that the guess provided to the BiCGSTAB is in fact the analytical solution to the problem, we do not observe significantly lower number of solver iterations.

Nevertheless, from a practical point of view, solve time is more important than the number of iterations required to obtain a numerical solution. The last column of the above figures clearly shows that the two direct solvers significantly outperform the iterative BiCGSTAB solver. On the other hand, employing BiCGSTAB with a ILUT preconditioner or providing guess to it, seems to have marginal effect on the required time.

== The three-dimensional Poisson problem ==

Next, we test the solvers on a three-dimensional Poisson problem. The governing problem is exactly the same as in the previous case, except that the domain is a $d=3$ dimensional unit sphere.

=== Solver performance ===

Below, we show the performance of the selected solvers in case of a two-dimensional Poisson problem. The first column shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of iterations a given solver required (direct solvers were assigned a value 1) and the third column shows the total time required to solve the system and obtain the numerical solution $\widehat u$.

[[File:convergence_3d_poisson_sin.png|1200px]]

[[File:convergence_3d_poisson_exp.png|1200px]]

In terms of accuracy, we again observe some unexpected behaviour when ILUT preconditioner is used. This time, the ILUT preconditioner was used to limit the solver iterations to 300 with fill factor 40, drop tolerance $10^{-5}$ and global tolerance $10^{-15}$. Although 100 more iterations were allowed than in the previous two-dimensional case, this clearly was not enough for the finest domain discretizations. Additionally, providing BiCGSTAB with the analytical guess yields even more confusing results as it yielded a sudden improve of the solution accuracy. This is, of course, wrong. It is important to note that we provide '''analytical solution''' as the guess. Therefore, as soon as the solver hits the maximum allowed iterations threshold it returns the best solution currently available - largely consisting out of the analytical guess provided.

In fact, the most reliable results are again yielded by the pure BiCGSTAB and the two direct solvers. However, the two direct solvers both ran out of the available 11.5 Gb RAM space available. The SparseLU was still able to obtain a solution with $10^5$ discretization nodes, while PardisoLU seems to be more efficient since it was still able to obtain a solution with $4\cdot 10^5$ discretization points.

From the time perspective, it is clear that the two direct solvers were again significantly faster than the iterative solver - that is, of course, until they ran out of RAM. We did not come across any RAM related
issues with the BiCGSTAB, even when the domain was discretized with more than a million nodes.

== The linear elasticity: Boussinesq's problem ==

As a final problem we solve the three-dimensional Boussinesq's problem, where a concentrated normal traction acts on an isotropic half-space <ref name="slaughter">Slaughter, William S. The linearized theory of elasticity. Springer Science & Business Media, 2012.</ref>. The same problem has also just recently been solved by employing $h$-adaptive solution procedure <ref name="slak">J. Slak and G. Kosec, Adaptive radial basis function-generated finite differences method for contact problems, International Journal for Numerical Methods in Engineering, 2019. DOI: 10.1002/nme.6067.</ref>, to demonstrate the performance of a hybrid mesh-free approximation method<ref name="hybrid">M. Jančič and G. Kosec, "A hybrid RBF-FD and WLS mesh-free strong-form approximation method," 2022 7th International Conference on Smart and Sustainable Technologies (SpliTech), 2022, pp. 1-6, doi: 10.23919/SpliTech55088.2022.9854278.</ref> also to test the performance of $hp$-adaptive solution procedures <ref name="hp">Jančič, Mitja, and Gregor Kosec. "Strong form mesh-free $ hp $-adaptive solution of linear elasticity problem." arXiv preprint arXiv:2210.07073 (2022).</ref>.

The problem has a closed form solution given in cylindrical coordinates $r$, $\theta$ and $z$ as

$$u_r = \frac{Pr}{4\pi \mu} \left(\frac{z}{R^3} - \frac{1-2\nu}{R(z+R)} \right), \qquad u_\theta = 0, \qquad u_z =
\frac{P}{4\pi \mu} \left(\frac{2(1-\nu)}{R} + \frac{z^2}{R^3} \right)$$

with stress tensor components

$$ \sigma_{rr} = \frac{P}{2\pi} \left( \frac{1-2\nu}{R(z+R)} - \frac{3r^2z}{R^5} \right), \qquad \sigma_{\theta\theta} =
\frac{P(1-2\nu)}{2\pi} \left( \frac{z}{R^3} - \frac{1}{R(z+R)} \right), \qquad \sigma_{zz} = -\frac{3Pz^3}{2 \pi R^5},$$

$$ \sigma_{rz} = -\frac{3Prz^2}{2 \pi R^5} , \quad \sigma_{r\theta} = 0 , \qquad\sigma_{\theta z} = 0$$

where $P$ is the magnitude of the concentrated force, $\nu$ is the Poisson's ratio, $\mu$ is the Lame parameter
and $R$ is the Eucledian distance to the origin. The solution has a singularity at the origin, which makes the problem ideal for treatment with adaptive procedures. Furthermore, the closed form solution also allows us to evaluate the accuracy of the numerical solution.

In our setup, we consider only a small part of the problem, i.e., $\varepsilon = 0.1$ away from the singularity, as schematically shown in figure below. From a numerical point of view, we solve the Navier-Cauchy Equation with Dirichlet boundary conditions described in the sketch, where the domain $\Omega$ is defined as a box, i.e., $\Omega = [-1, -\varepsilon] \times [-1, -\varepsilon] \times [-1, -\varepsilon]$.

[[File:boussinesq_sketch.png|400px]]

Example solution with approximately $2\cdot 10^5$ nodes is shown below.

[[File:solution_boussinesq.png|400px]]

=== Solver performance ===

Below, we show the performance of the selected solvers in case of a three-dimensional Boussinesq's problem. The first column shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of iterations a given solver required (direct solvers were assigned a value 1) and the third column shows the total time required to solve the system and obtain the numerical solution $\widehat u$.

[[File:convergence_3d_boussinesq.png|1200px]]

Similarly to previous studies on Poisson problems, we observe that in terms of accuracy of the numerical solution, all chosen solvers are in good agreement. However, this time we observe some desirable effects of the ILUT and guess on the performance of BiCGSTAB. Using ILUT preconditioner (and) guess reduces the number of iterations required - however, this does not notably improve the times required to obtain a numerical solution. Therefore, considering the fact that a badly-thought-out preconditioner orguess definition can notably affect the solver's performance, one is probably safer using the general purpose BiCGSTAB without the guess and without the preconditioner.

The two direct solvers were again able to solve the linear sparse system notably faster than the iterative BiCGSTAB. Unfortunately, using direct solvers we again ran into RAM space related issues.

= Further literature =

To read more about the solvers, we propose:
* [https://users.fmf.uni-lj.si/plestenjak/Vaje/INMLA/Predavanja/Skripta.pdf Slovenian script by Bor Plestenjak]
* Matrix Algorithms: Volume II: Eigensystems<ref name="gilbert1">Stewart, Gilbert W. Matrix Algorithms: Volume II: Eigensystems. Society for Industrial and Applied Mathematics, 2001.</ref>
* Matrix algorithms: volume 1: basic decompositions<ref name="gilbert2">Stewart, Gilbert W. Matrix algorithms: volume 1: basic decompositions. Society for Industrial and Applied Mathematics, 1998.</ref>

=References=
<references/>

Solving sparse systems

2023-02-14T12:56:57Z

Mitjajancic: Added some observations.

There are many methods available for solving sparse systems. We compare some of them here.
<figure id="fig:matrix1">
[[File:matrix.png|300px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Matrix of the discretized PDE. </caption>]]
</figure>
Mathematica has the following methods available (https://reference.wolfram.com/language/ref/LinearSolve.html#DetailsAndOptions)
* direct: banded, cholesky, multifrontal (direct sparse LU)
* iterative: Krylov

Matlab has the following methods:
* direct: https://www.mathworks.com/help/matlab/ref/mldivide.html#bt42omx_head
* iterative: https://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html#brzoiix, including bicgstab, gmres

Eigen has the following methods: (https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
* direct: sparse LU
* iterative: bicgstab, cg

Solving a simple sparse system $A x = b$ for steady space of heat equation in 1d with $n$ nodes, results in a matrix shown in <xr id="fig:matrix1"/>.

The following timings of solvers are given in seconds:
{| class="wikitable"
|-
! $n = 10^6$
! Matlab
! Mathematica
! Eigen
|-
! Banded
| 0.16
| 0.28
| 0.04
|-
! SparseLU
| /
| 1.73
| 0.82
|-
|}

Incomplete LU preconditioner was used for BICGStab.
Without the preconditioner BICGStab does not converge.

==Parallel execution ==
BICGStab can be run in parallel, as explain in the general parallelism: https://eigen.tuxfamily.org/dox/TopicMultiThreading.html, and specifically

'''"When using sparse matrices, best performance is achied for a row-major sparse matrix format. Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled".'''

Eigen uses number of threads specified my OopenMP, unless <code>Eigen::setNbThreads(n);</code> was called.
Minimal working example:

<figure id="fig:eigen_par">
[[File:eigen_parallel_2.png|1200px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Memory and CPU usage. C stands for construction of the system, L stands for the calculation of ILUT preconditioner and S for BICGStab iteration.</caption>]]
</figure>

<syntaxhighlight lang="c++" line>
#include <iostream>
#include <vector>
#include "Eigen/Sparse"
#include "Eigen/IterativeLinearSolvers"

using namespace std;
using namespace Eigen;

int main(int argc, char* argv[]) {

assert(argc == 2 && "Second argument is size of the system.");
stringstream ss(argv[1]);
int n;
ss >> n;
cout << "n = " << n << endl;

// Fill the matrix
VectorXd b = VectorXd::Ones(n) / n / n;
b(0) = b(n-1) = 1;
SparseMatrix<double, RowMajor> A(n, n);
A.reserve(vector<int>(n, 3)); // 3 per row
for (int i = 0; i < n-1; ++i) {
A.insert(i, i) = -2;
A.insert(i, i+1) = 1;
A.insert(i+1, i) = 1;
}
A.coeffRef(0, 0) = 1;
A.coeffRef(0, 1) = 0;
A.coeffRef(n-1, n-2) = 0;
A.coeffRef(n-1, n-1) = 1;

// Solve the system
BiCGSTAB<SparseMatrix<double, RowMajor>, IncompleteLUT<double>> solver;
solver.setTolerance(1e-10);
solver.setMaxIterations(1000);
solver.compute(A);
VectorXd x = solver.solve(b);
cout << "#iterations: " << solver.iterations() << endl;
cout << "estimated error: " << solver.error() << endl;
cout << "sol: " << x.head(6).transpose() << endl;

return 0;
}
</syntaxhighlight>

was compiled with <code>g++ -o parallel_solve -O3 -fopenmp solver_test_parallel.cpp</code>.
<xr id="fig:eigen_par"/> was produced when the program above was run as <code>./parallel_solve 10000000</code>.

==ILUT factorization==
This is a method to compute a general algebraic preconditioner which has the ability to control (to some degree) the memory and time complexity of the solution.
It is described in a paper by

Saad, Yousef. "ILUT: A dual threshold incomplete LU factorization." Numerical linear algebra with applications 1.4 (1994): 387-402.

It has two parameters, tolerance $\tau$ and fill factor $f$. Elements in $L$ and $U$ factors that would be below relative tolerance $\tau$ (say, to 2-norm of this row) are dropped.
Then, only $f$ elements per row in $L$ and $U$ are allowed and the maximum $f$ (if they exist) are taken. The diagonal elements are always kept. This means that the number of non-zero
elements in the rows of the preconditioner can not exceed $2f+1$.

Greater parameter $f$ means slover but better $ILUT$ factorization. If $\tau$ is small, having a very large $f$ means getting a complete $LU$ factorization. Typical values are $5, 10, 20, 50$.
Lower $\tau$ values mean that move small elements are kept, resulting in a more precise, but longer computation. With a large $f$, the factorization approaches the direct $LU$ factorization as $\tau \to 0$. Typical values are from 1e-2
to 1e-6.

Whether the method converges or diverges for given parameters is very sudden and should be determined for each apllication separately. An example of
such a convergence graph for the diffusion equation is presented in figures below:

[[File:bicgstab_err.png|400px]][[File:bicgstab_conv.png|400px]][[File:bicgstab_iter.png|400px]]

We can see that low BiCGSTAB error estimate coincides with the range of actual small error compared to the analytical solution as well as with the range of lesser number of iterations.

==Pardiso library==
[[Pardiso]]

==Herzian contact ==
The matrix of the [[Hertzian contact]] problem is shown below:
[[File:matrix_hertzian.png|500px]]

The condition number iz arther large and grows superlinearly with size in the example below.
[[File:cond.png|600px]]

= Discussion based on practical observations =

In this section we analyse different solvers, all available as part of
the [https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html Eigen library] supported by
the [https://e6.ijs.si/medusa/ Medusa]. Namely we analyse the performance
of two direct solvers

* [https://eigen.tuxfamily.org/dox/classEigen_1_1SparseLU.html SparseLU]
* [https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html SparseQR]

three combinations of general purpose iterative
solver [https://eigen.tuxfamily.org/dox/classEigen_1_1BiCGSTAB.html BiCGSTAB]

* BiCGSTAB
* BiCGSTAB with ILUT preconditioner
* BiCGSTAB with ILUT preconditioner and exploited '''[https://eigen.tuxfamily.org/dox/classEigen_1_1SolveWithGuess.html solveWithGuess()]''' functionality

and a wrapper to external direct solver

* [https://eigen.tuxfamily.org/dox/classEigen_1_1PardisoLU.html PardisoLU].

The performance tests are focused on the accuracy of the numerical solution, number of iterations and time required to solve the final sparse system. We also comment on the amount of RAM required by a specific solver. If required, the [https://gitlab.com/e62Lab/eigen_solvers_analysis gitlab repository] with the implementation is also provided.

<syntaxhighlight lang="text">
NOTE: Although the SparseQR has been included in our tests, significantly larger solve times compared to other solvers have been observed. For that reason, the SparseQR is discarded from further analyses even though it was able to obtain a solution of comparable accuracy.
</syntaxhighlight>

== The two-dimensional Poisson problem ==

First, we solve the Poisson problem
$$\nabla ^2 u(\boldsymbol x) = f(\boldsymbol x),$$

for $\boldsymbol x \in \Omega $ where the domain $\Omega$ is a unit $d=2$ dimensional circle with radius $R=1$.
We use mixed boundary conditions, i.e., Dirichlet boundary for all $x < 0$ and Neumann boundary for all $x \geq 0$. Additionally, ghost or fictitious nodes were added to the Neumann boundary.

We perform analyses on two scenarios, i.e., two definitions of $f(\boldsymbol x)$:

* $f(\boldsymbol x) = -d\pi ^2 \prod_{i = 1}^d \sin(\pi x_i)$ yielding
** $f_d(\boldsymbol x) = \prod_{i = 1}^d \sin(\pi x_i)$ on the Dirichlet boundary and
** $ \boldsymbol f_n(\boldsymbol x) = \pi \sum_{i = 1}^d \cos(\pi x_i) \prod_{j=1, j\neq i}^d \sin (\pi x_j)
\boldsymbol e_i$ with unit vectors $\boldsymbol e_i$ on the Neumann boundary.
* $f(\boldsymbol x) = 2a e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}(2a\left \| \boldsymbol x -
\boldsymbol x_s \right \| - d)$ yielding
** $f_d(\boldsymbol x) = e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$ yielding on Dirichlet boundary and
** $\boldsymbol f_n (\boldsymbol x) = -2a(\boldsymbol x - \boldsymbol x_s)e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$ on Neumann boundary.

for width $a=10^2$ of the exponential source positioned at $\boldsymbol x_s=\boldsymbol 0$.

Example solutions to both scenarios are shown in figure below.

[[File:solution_poisson.png|800px]]

=== Solver performance ===

Below, we show the performance of the selected solvers in case of a two-dimensional Poisson problem. The first column shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of iterations a given solver required (direct solvers were assigned a value 1) and the third column shows the total time required to solve the system and obtain the numerical solution $\widehat u$.

<syntaxhighlight lang="text">
Computational times: To obtain the computational times the code was run on a single CPU with assured CPU affinity using the ''taskset'' command.
</syntaxhighlight>

[[File:convergence_poisson_sin.png|1200px]]

[[File:convergence_poisson_exp.png|1200px]]

In terms of accuracy, all solvers under consideration perform equally good, with an exception of BiCGSTAB with ILUT preconditioner, who failed to obtain a solution of sufficient accuracy at the finest domain discretization with approximately $2\cdot 10^6$ nodes. The reason is most likely in the ILUT preconditioner, set to limit solver iterations to 200 with fill factor 40, drop tolerance $10^{-5}$ and global tolerance $10^{-15}$. A different ILUT setup would most probably yield a more desirable behaviour, but it nevertheless shows that when using a ILUT preconditioner one must be careful before making conclusions. It is also important to note that the SparseLU ran out of the available 11.5 Gb RAM and, therefore, failed to obtain a numerical solution. Other solvers had no issues with the available RAM space.

It is perhaps more interesting to study the number of iterations required by the solvers. Obviously, direct solvers require a "single iteration". Iterative solvers, on the other hand, require multiple iterations. Generally speaking, we observe that the number of solver iterations increases with the number of discretization points. However, it is interesting to see that the pure BiCGSTAB shows the most stable and predictable behaviour. Despite the fact that the guess provided to the BiCGSTAB is in fact the analytical solution to the problem, we do not observe significantly lower number of solver iterations.

Nevertheless, from a practical point of view, solve time is more important than the number of iterations required to obtain a numerical solution. The last column of the above figures clearly shows that the two direct solvers significantly outperform the iterative BiCGSTAB solver. On the other hand, employing BiCGSTAB with a ILUT preconditioner or providing guess to it, seems to have marginal effect on the required time.

== The three-dimensional Poisson problem ==

Next, we test the solvers on a three-dimensional Poisson problem. The governing problem is exactly the same as in the previous case, except that the domain is a $d=3$ dimensional unit sphere.

=== Solver performance ===

Below, we show the performance of the selected solvers in case of a two-dimensional Poisson problem. The first column shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of iterations a given solver required (direct solvers were assigned a value 1) and the third column shows the total time required to solve the system and obtain the numerical solution $\widehat u$.

[[File:convergence_3d_poisson_sin.png|1200px]]

[[File:convergence_3d_poisson_exp.png|1200px]]

In terms of accuracy, we again observe some unexpected behaviour when ILUT preconditioner is used. This time, the ILUT preconditioner was used to limit the solver iterations to 300 with fill factor 40, drop tolerance $10^{-5}$ and global tolerance $10^{-15}$. Although 100 more iterations were allowed than in the previous two-dimensional case, this clearly was not enough for the finest domain discretizations. Additionally, providing BiCGSTAB with the analytical guess yields even more confusing results as it yielded a sudden improve of the solution accuracy. This is, of course, wrong. It is important to note that we provide '''analytical solution''' as the guess. Therefore, as soon as the solver hits the maximum allowed iterations threshold it returns the best solution currently available - largely consisting out of the analytical guess provided.

In fact, the most reliable results are again yielded by the pure BiCGSTAB and the two direct solvers. However, the two direct solvers both ran out of the available 11.5 Gb RAM space available. The SparseLU was still able to obtain a solution with $10^5$ discretization nodes, while PardisoLU seems to be more efficient since it was still able to obtain a solution with $4\cdot 10^5$ discretization points.

From the time perspective, it is clear that the two direct solvers were again significantly faster than the iterative solver - that is, of course, until they ran out of RAM. We did not come across any RAM related
issues with the BiCGSTAB, even when the domain was discretized with more than a million nodes.

== The linear elasticity: Boussinesq's problem ==

As a final problem we solve the three-dimensional Boussinesq's problem, where a concentrated normal traction acts on an isotropic half-space <ref name="slaughter">Slaughter, William S. The linearized theory of elasticity. Springer Science & Business Media, 2012.</ref>. The same problem has also just recently been solved by employing $h$-adaptive solution procedure <ref name="slak">J. Slak and G. Kosec, Adaptive radial basis function-generated finite differences method for contact problems, International Journal for Numerical Methods in Engineering, 2019. DOI: 10.1002/nme.6067.</ref>, to demonstrate the performance of a hybrid mesh-free approximation method<ref name="hybrid">M. Jančič and G. Kosec, "A hybrid RBF-FD and WLS mesh-free strong-form approximation method," 2022 7th International Conference on Smart and Sustainable Technologies (SpliTech), 2022, pp. 1-6, doi: 10.23919/SpliTech55088.2022.9854278.</ref> also to test the performance of $hp$-adaptive solution procedures <ref name="hp">Jančič, Mitja, and Gregor Kosec. "Strong form mesh-free $ hp $-adaptive solution of linear elasticity problem." arXiv preprint arXiv:2210.07073 (2022).</ref>.

The problem has a closed form solution given in cylindrical coordinates $r$, $\theta$ and $z$ as

$$u_r = \frac{Pr}{4\pi \mu} \left(\frac{z}{R^3} - \frac{1-2\nu}{R(z+R)} \right), \qquad u_\theta = 0, \qquad u_z =
\frac{P}{4\pi \mu} \left(\frac{2(1-\nu)}{R} + \frac{z^2}{R^3} \right)$$

with stress tensor components

$$ \sigma_{rr} = \frac{P}{2\pi} \left( \frac{1-2\nu}{R(z+R)} - \frac{3r^2z}{R^5} \right), \qquad \sigma_{\theta\theta} =
\frac{P(1-2\nu)}{2\pi} \left( \frac{z}{R^3} - \frac{1}{R(z+R)} \right), \qquad \sigma_{zz} = -\frac{3Pz^3}{2 \pi R^5},$$

$$ \sigma_{rz} = -\frac{3Prz^2}{2 \pi R^5} , \quad \sigma_{r\theta} = 0 , \qquad\sigma_{\theta z} = 0$$

where $P$ is the magnitude of the concentrated force, $\nu$ is the Poisson's ratio, $\mu$ is the Lame parameter
and $R$ is the Eucledian distance to the origin. The solution has a singularity at the origin, which makes the problem ideal for treatment with adaptive procedures. Furthermore, the closed form solution also allows us to evaluate the accuracy of the numerical solution.

In our setup, we consider only a small part of the problem, i.e., $\varepsilon = 0.1$ away from the singularity, as schematically shown in figure below. From a numerical point of view, we solve the Navier-Cauchy Equation with Dirichlet boundary conditions described in the sketch, where the domain $\Omega$ is defined as a box, i.e., $\Omega = [-1, -\varepsilon] \times [-1, -\varepsilon] \times [-1, -\varepsilon]$.

[[File:boussinesq_sketch.png|400px]]

Example solution with approximately $2\cdot 10^5$ nodes is shown below.

[[File:solution_boussinesq.png|400px]]

=== Solver performance ===

Below, we show the performance of the selected solvers in case of a three-dimensional Boussinesq's problem. The first column shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of iterations a given solver required (direct solvers were assigned a value 1) and the third column shows the total time required to solve the system and obtain the numerical solution $\widehat u$.

[[File:convergence_3d_boussinesq.png|1200px]]

Similarly to previous studies on Poisson problems, we observe that in terms of accuracy of the numerical solution, all chosen solvers are in good agreement. However, this time we observe some desirable effects of the ILUT and guess on the performance of BiCGSTAB. Using ILUT preconditioner (and) guess reduces the number of iterations required - however, this does not notably improve the times required to obtain a numerical solution. Therefore, considering the fact that a badly-thought-out preconditioner orguess definition can notably affect the solver's performance, one is probably safer using the general purpose BiCGSTAB without the guess and without the preconditioner.

The two direct solvers were again able to solve the linear sparse system notably faster than the iterative BiCGSTAB. Unfortunately, using direct solvers we again ran into RAM space related issues.

= Further literature =

To read more about the solvers, we propose:
* [https://users.fmf.uni-lj.si/plestenjak/Vaje/INMLA/Predavanja/Skripta.pdf Slovenian script by Bor Plestenjak]
* Matrix Algorithms: Volume II: Eigensystems<ref name="gilbert1">Stewart, Gilbert W. Matrix Algorithms: Volume II: Eigensystems. Society for Industrial and Applied Mathematics, 2001.</ref>
* Matrix algorithms: volume 1: basic decompositions<ref name="gilbert2">Stewart, Gilbert W. Matrix algorithms: volume 1: basic decompositions. Society for Industrial and Applied Mathematics, 1998.</ref>

=References=
<references/>

File:Convergence 3d boussinesq.png

2023-02-14T12:36:47Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Solution boussinesq.png

2023-02-14T12:35:53Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Boussinesq sketch.png

2023-02-14T12:35:53Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Convergence 3d poisson sin.png

2023-02-14T12:32:30Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Convergence 3d poisson exp.png

2023-02-14T12:32:30Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Convergence poisson exp.png

2023-02-14T12:30:17Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Convergence poisson sin.png

2023-02-14T12:30:17Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Solution poisson.png

2023-02-14T12:28:19Z

Mitjajancic: Mitjajancic uploaded a new version of File:Solution poisson.png

File uploaded with MsUpload

File:Solution poisson.png

2023-02-14T12:27:47Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

Solving sparse systems

2023-02-14T12:27:34Z

Mitjajancic:

There are many methods available for solving sparse systems. We compare some of them here.
<figure id="fig:matrix1">
[[File:matrix.png|300px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Matrix of the discretized PDE. </caption>]]
</figure>
Mathematica has the following methods available (https://reference.wolfram.com/language/ref/LinearSolve.html#DetailsAndOptions)
* direct: banded, cholesky, multifrontal (direct sparse LU)
* iterative: Krylov

Matlab has the following methods:
* direct: https://www.mathworks.com/help/matlab/ref/mldivide.html#bt42omx_head
* iterative: https://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html#brzoiix, including bicgstab, gmres

Eigen has the following methods: (https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
* direct: sparse LU
* iterative: bicgstab, cg

Solving a simple sparse system $A x = b$ for steady space of heat equation in 1d with $n$ nodes, results in a matrix shown in <xr id="fig:matrix1"/>.

The following timings of solvers are given in seconds:
{| class="wikitable"
|-
! $n = 10^6$
! Matlab
! Mathematica
! Eigen
|-
! Banded
| 0.16
| 0.28
| 0.04
|-
! SparseLU
| /
| 1.73
| 0.82
|-
|}

Incomplete LU preconditioner was used for BICGStab.
Without the preconditioner BICGStab does not converge.

==Parallel execution ==
BICGStab can be run in parallel, as explain in the general parallelism: https://eigen.tuxfamily.org/dox/TopicMultiThreading.html, and specifically

'''"When using sparse matrices, best performance is achied for a row-major sparse matrix format. Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled".'''

Eigen uses number of threads specified my OopenMP, unless <code>Eigen::setNbThreads(n);</code> was called.
Minimal working example:

<figure id="fig:eigen_par">
[[File:eigen_parallel_2.png|1200px|thumb|upright=2|alt=Matrix of the discretized PDE.|<caption>Memory and CPU usage. C stands for construction of the system, L stands for the calculation of ILUT preconditioner and S for BICGStab iteration.</caption>]]
</figure>

<syntaxhighlight lang="c++" line>
#include <iostream>
#include <vector>
#include "Eigen/Sparse"
#include "Eigen/IterativeLinearSolvers"

using namespace std;
using namespace Eigen;

int main(int argc, char* argv[]) {

assert(argc == 2 && "Second argument is size of the system.");
stringstream ss(argv[1]);
int n;
ss >> n;
cout << "n = " << n << endl;

// Fill the matrix
VectorXd b = VectorXd::Ones(n) / n / n;
b(0) = b(n-1) = 1;
SparseMatrix<double, RowMajor> A(n, n);
A.reserve(vector<int>(n, 3)); // 3 per row
for (int i = 0; i < n-1; ++i) {
A.insert(i, i) = -2;
A.insert(i, i+1) = 1;
A.insert(i+1, i) = 1;
}
A.coeffRef(0, 0) = 1;
A.coeffRef(0, 1) = 0;
A.coeffRef(n-1, n-2) = 0;
A.coeffRef(n-1, n-1) = 1;

// Solve the system
BiCGSTAB<SparseMatrix<double, RowMajor>, IncompleteLUT<double>> solver;
solver.setTolerance(1e-10);
solver.setMaxIterations(1000);
solver.compute(A);
VectorXd x = solver.solve(b);
cout << "#iterations: " << solver.iterations() << endl;
cout << "estimated error: " << solver.error() << endl;
cout << "sol: " << x.head(6).transpose() << endl;

return 0;
}
</syntaxhighlight>

was compiled with <code>g++ -o parallel_solve -O3 -fopenmp solver_test_parallel.cpp</code>.
<xr id="fig:eigen_par"/> was produced when the program above was run as <code>./parallel_solve 10000000</code>.

==ILUT factorization==
This is a method to compute a general algebraic preconditioner which has the ability to control (to some degree) the memory and time complexity of the solution.
It is described in a paper by

Saad, Yousef. "ILUT: A dual threshold incomplete LU factorization." Numerical linear algebra with applications 1.4 (1994): 387-402.

It has two parameters, tolerance $\tau$ and fill factor $f$. Elements in $L$ and $U$ factors that would be below relative tolerance $\tau$ (say, to 2-norm of this row) are dropped.
Then, only $f$ elements per row in $L$ and $U$ are allowed and the maximum $f$ (if they exist) are taken. The diagonal elements are always kept. This means that the number of non-zero
elements in the rows of the preconditioner can not exceed $2f+1$.

Greater parameter $f$ means slover but better $ILUT$ factorization. If $\tau$ is small, having a very large $f$ means getting a complete $LU$ factorization. Typical values are $5, 10, 20, 50$.
Lower $\tau$ values mean that move small elements are kept, resulting in a more precise, but longer computation. With a large $f$, the factorization approaches the direct $LU$ factorization as $\tau \to 0$. Typical values are from 1e-2
to 1e-6.

Whether the method converges or diverges for given parameters is very sudden and should be determined for each apllication separately. An example of
such a convergence graph for the diffusion equation is presented in figures below:

[[File:bicgstab_err.png|400px]][[File:bicgstab_conv.png|400px]][[File:bicgstab_iter.png|400px]]

We can see that low BiCGSTAB error estimate coincides with the range of actual small error compared to the analytical solution as well as with the range of lesser number of iterations.

==Pardiso library==
[[Pardiso]]

==Herzian contact ==
The matrix of the [[Hertzian contact]] problem is shown below:
[[File:matrix_hertzian.png|500px]]

The condition number iz arther large and grows superlinearly with size in the example below.
[[File:cond.png|600px]]

= Discussion based on practical observations =

In this section we analyse different solvers, all available as part of
the [https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html Eigen library] supported by
the [https://e6.ijs.si/medusa/ Medusa]. Namely we analyse the performance
of two direct solvers

* [https://eigen.tuxfamily.org/dox/classEigen_1_1SparseLU.html SparseLU]
* [https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html SparseQR]

three combinations of general purpose iterative
solver [https://eigen.tuxfamily.org/dox/classEigen_1_1BiCGSTAB.html BiCGSTAB]

* BiCGSTAB
* BiCGSTAB with ILUT preconditioner
* BiCGSTAB with ILUT preconditioner and exploited '''[https://eigen.tuxfamily.org/dox/classEigen_1_1SolveWithGuess.html solveWithGuess()]''' functionality

and a wrapper to external direct solver

* [https://eigen.tuxfamily.org/dox/classEigen_1_1PardisoLU.html PardisoLU].

The performance tests are focused on the accuracy of the numerical solution, number of iterations and time required to solve the final sparse system. We also comment on the amount of RAM required by a specific solver.

<syntaxhighlight lang="text">
NOTE: Although the SparseQR has been included in our tests, significantly larger solve times compared to other solvers have been observed. For that reason, the SparseQR is discarded from further analyses even though it was able to obtain a solution of comparable accuracy.
</syntaxhighlight>

== The two-dimensional Poisson problem ==

First, we solve the Poisson problem
$$\nabla ^2 u(\boldsymbol x) = f(\boldsymbol x),$$

for $\boldsymbol x \in \Omega $ where the domain $\Omega$ is a unit $d=2$ dimensional circle with radius $R=1$.
We use mixed boundary conditions, i.e., Dirichlet boundary for all $x < 0$ and Neumann boundary for all $x \geq 0$. Additionally, ghost or fictitious nodes were added to the Neumann boundary.

We perform analyses on two scenarios, i.e., two definitions of $f(\boldsymbol x)$:

* $f(\boldsymbol x) = -d\pi ^2 \prod_{i = 1}^d \sin(\pi x_i)$ yielding
** $f_d(\boldsymbol x) = \prod_{i = 1}^d \sin(\pi x_i)$ on the Dirichlet boundary and
** $ \boldsymbol f_n(\boldsymbol x) = \pi \sum_{i = 1}^d \cos(\pi x_i) \prod_{j=1, j\neq i}^d \sin (\pi x_j)
\boldsymbol e_i$ with unit vectors $\boldsymbol e_i$ on the Neumann boundary.
* $f(\boldsymbol x) = 2a e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}(2a\left \| \boldsymbol x -
\boldsymbol x_s \right \| - d)$ yielding
** $f_d(\boldsymbol x) = e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$ yielding on Dirichlet boundary and
** $\boldsymbol f_n (\boldsymbol x) = -2a(\boldsymbol x - \boldsymbol x_s)e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$ on Neumann boundary.

for width $a=10^2$ of the exponential source positioned at $\boldsymbol x_s=\boldsymbol 0$.

Example solutions to both scenarios are shown in figure below.

[[File:solution_poisson.jpg]]

### Solver performance

Below, we show the performance of the selected solvers in case of a two-dimensional Poisson problem. The first column
shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of
iterations a given solver required (direct solvers were assigned a value 1) and the third column shows
the total time required to solve the system and obtain the numerical solution $'\widehat u'$.

> **Computational times:**
>
>To obtain the computational times the code was run on a single CPU with assured CPU affinity using the *taskset*
> command.

![](../results/convergence_poisson_sin.png)
![](../results/convergence_poisson_exp.png)

In terms of accuracy, all solvers under consideration perform equally good, with an exception of BiCGSTAB with ILUT
preconditioner,
who failed to obtain a solution of sufficient accuracy at the finest domain discretization with approximately $'2\cdot
10^6'$ nodes. The reason is most likely
in the ILUT preconditioner, set to limit solver iterations to 200 with fill factor 40, drop tolerance $'10^{-5}'$ and
global tolerance $'10^{-15}'$. A different ILUT setup would most probably yield a more desirable behaviour, but it
nevertheless shows that when using a ILUT preconditioner one must be careful before making conclusions. It is also
important to note that the SparseLU ran out of the available 11.5 Gb RAM and, therefore, failed to obtain a numerical
solution. Other solvers had no issues with the available RAM space.

It is perhaps more interesting to study the number of iterations required by the solvers. Obviously, direct solvers
require a "single iteration". Iterative solvers, on the other hand, require multiple iterations. Generally speaking, we
observe that the number of solver iterations increases with the number of discretization points. However, it is
interesting to see that the pure BiCGSTAB shows the most stable and predictable behaviour. Despite the fact that the
guess provided to the BiCGSTAB is in fact the analytical solution to the problem, we do not observe significantly lower
number of solver iterations.

Nevertheless, from a practical point of view, solve time is more important than the number of iterations required to
obtain a numerical solution. The last column of the above figures clearly shows that the two direct solvers
significantly outperform the iterative BiCGSTAB solver. On the other hand, employing BiCGSTAB with a ILUT preconditioner
or providing guess to it, seems to have marginal effect on the required time.

## The three-dimensional Poisson problem

Next, we test the solvers on a three-dimensional Poisson problem. The governing problem is exactly the same as in the
previous case, except that the domain is a $'d=3'$ dimensional unit sphere.

### Solver performance

Below, we show the performance of the selected solvers in case of a two-dimensional Poisson problem. The first column
shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of
iterations a given solver required (direct solvers were assigned a value 1) and the third column shows
the total time required to solve the system and obtain the numerical solution $'\widehat u'$.

![](../results/convergence_3d_poisson_sin.png)
![](../results/convergence_3d_poisson_exp.png)

In terms of accuracy, we again observe some unexpected behaviour when ILUT preconditioner is used. This time, the ILUT
preconditioner was used to limit the solver iterations to 300 with fill factor 40, drop tolerance $'10^{-5}'$ and global
tolerance $'10^{-15}'$. Although 100 more iterations were allowed than in the previous two-dimensional case, this
clearly was not enough for the finest domain discretizations. Additionally, providing BiCGSTAB with the analytical guess
yields even more confusing results as it yielded a sudden improve of the solution accuracy. This is, of
course, wrong. It is important to note that we provide **analytical solution** as the guess. Therefore, as soon as the
solver hits the
maximum allowed iterations threshold it returns the best solution currently available - largely consisting out of the
analytical guess provided.

In fact, the most reliable results are again yielded by the pure BiCGSTAB and the two direct solvers. However, the two
direct solvers both ran out of the available 11.5
Gb RAM space available. The SparseLU was still able to obtain a solution with $'10^5'$ discretization nodes, while
PardisoLU seems to be more efficient since it was still able to obtain a solution with $'4\cdot 10^5'$ discretization
points.

From the time perspective, it is clear that the two direct solvers were again significantly
faster than the iterative solver - that is, of course, until they ran out of RAM. We did not come across any RAM related
issues with the BiCGSTAB, even when the domain was discretized with more than a million nodes.

## The linear elasticity: Boussinesq's problem

As a final problem we solve the three-dimensional Boussinesq's problem, where a concentrated normal traction acts on an
isotropic half-space~\cite{Slaughter_2002}. The same problem has also just recently been solved by employing $'h'
$-adaptive
solution procedure~\cite{slak2019adaptive}, to demonstrate the performance of a hybrid mesh-free approximation
method~\cite{jancichybrid} also to test the performance of $'hp'$-adaptive solution procedures~\cite{hp}.

The problem has a closed form solution given in cylindrical coordinates $'r'$, $'\theta'$ and $'z'$ as

$'u_r = \frac{Pr}{4\pi \mu} \left(\frac{z}{R^3} - \frac{1-2\nu}{R(z+R)} \right), \qquad u_\theta = 0, \qquad u_z =
\frac{P}{4\pi \mu} \left(\frac{2(1-\nu)}{R} + \frac{z^2}{R^3} \right)'$

with stress tensor components

$' \sigma_{rr} = \frac{P}{2\pi} \left( \frac{1-2\nu}{R(z+R)} - \frac{3r^2z}{R^5} \right), \qquad \sigma_{\theta\theta} =
\frac{P(1-2\nu)}{2\pi} \left( \frac{z}{R^3} - \frac{1}{R(z+R)} \right), \qquad \sigma_{zz} = -\frac{3Pz^3}{2 \pi R^5}'$

$' \sigma_{rz} = -\frac{3Prz^2}{2 \pi R^5} , \quad \sigma_{r\theta} = 0 , \qquad\sigma_{\theta z} = 0'$

where $'P'$ is the magnitude of the concentrated force, $'\nu'$ is the Poisson's ratio, $'\mu'$ is the Lame parameter
and $R$ is the Eucledian distance to the origin. The solution has a singularity at the origin, which makes the problem
ideal for treatment with adaptive procedures. Furthermore, the closed form solution also allows us to evaluate the
accuracy of the numerical solution.

In our setup, we consider only a small part of the problem, i.e., $'\varepsilon = 0.1'$ away from the singularity, as
schematically shown in figure below. From a numerical point of view, we solve the Navier-Cauchy Equation with Dirichlet
boundary conditions described in the sketch, where the domain $'\Omega'$ is defined as a box, i.e., $'\Omega
= [-1, -\varepsilon] \times [-1, -\varepsilon] \times [-1, -\varepsilon]'$.

![](../results/boussinesq_sketch.png)

Example solution with approximately $'2\cdot 10^5'$ nodes is shown below.

![](../results/solution_boussinesq.png)

### Solver performance

Below, we show the performance of the selected solvers in case of a three-dimensional Boussinesq's problem. The first
column
shows the accuracy of the numerical solution in terms of the infinity norm, the second column shows the number of
iterations a given solver required (direct solvers were assigned a value 1) and the third column shows
the total time required to solve the system and obtain the numerical solution $'\widehat u'$.

![](../results/convergence_3d_boussinesq.png)

Similarly to previous studies on Poisson problems, we observe that in terms of accuracy of the numerical solution, all
chosen solvers are in good agreement.
However, this time we observe some desirable effects of the ILUT and guess on the performance of BiCGSTAB. Using ILUT
preconditioner (and) guess reduces the number of iterations required - however, this does not notably improve the times
required to obtain a numerical solution. Therefore, considering the fact that a badly-thought-out preconditioner or
guess definition can notably affect the solver's performance, one is probably safer using the general purpose
BiCGSTAB without the guess and without the preconditioner.

The two direct solvers were again able to solve the linear sparse system notably faster than the iterative BiCGSTAB.
Unfortunately, using direct solvers we again ran into RAM space related issues.

# More on solvers

- According to our observations, Pardiso is the fastest and most reliable, given the available RAM is large enough.
- https://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html

# Further literature:

Script: https://users.fmf.uni-lj.si/plestenjak/Vaje/INMLA/Predavanja/Skripta.pdf

Stewart, Gilbert W. Matrix Algorithms: Volume II: Eigensystems. Society for Industrial and Applied Mathematics, 2001.

Stewart, Gilbert W. Matrix algorithms: volume 1: basic decompositions. Society for Industrial and Applied Mathematics,
1998.

File:Fwo abaqus.png

2022-10-13T13:30:41Z

Mitjajancic: Mitjajancic uploaded a new version of File:Fwo abaqus.png

File uploaded with MsUpload

File:Error indicator convergence.png

2022-10-13T13:30:07Z

Mitjajancic: Mitjajancic uploaded a new version of File:Error indicator convergence.png

File uploaded with MsUpload

File:Refinement demonstration 2d.png

2022-10-13T13:28:48Z

Mitjajancic: Mitjajancic uploaded a new version of File:Refinement demonstration 2d.png

File uploaded with MsUpload

File:Refinement workflow.png

2022-10-13T13:26:41Z

Mitjajancic: Mitjajancic uploaded a new version of File:Refinement workflow.png

File uploaded with MsUpload

HP-adaptivity

2022-09-25T07:51:08Z

Mitjajancic:

Adaptive solution procedures are essential in problems where the accuracy of the numerical solution varies spatially and are currently subject of intensive studies. Two conceptually different adaptive approaches have been proposed, namely $p$-adaptivity or $h$-, $r$-adaptivity. In $p$-adaptivity, the accuracy of the numerical solution is varied by changing the order of approximation, while in $h$- and $r$-adaptivity the resolution of the spatial discretization is adjusted for the same purpose. In the $h$-adaptive approach, nodes are added or removed from the domain as needed, while in the $r$-adaptive approach, the total number of nodes remains constant - the nodes are only repositioned with respect to the desired accuracy. Ultimately, $h$- and $p$-adaptivities can be combined to form the so-called $hp$-adaptivity, where the accuracy of the solution is controlled with the order of the method and the resolution of the spatial discretization.

__FORCETOC__

== Basic concept ==

The proposed $hp$-adaptive solution procedure follows the well-established paradigm based on an iterative loop, where each iteration step consists of four modules:
* '''Solve''' - A numerical solution $\widehat{u}$ is obtained.
* '''Estimate''' - An estimate of the spatial accuracy of the numerical solution is calculated using error indicators.
* '''Mark''' - Depending on the error indicator values $\eta _i$, a marking strategy is used to mark the computational nodes for (de)refinement.
* '''Refine''' - Refinement strategy is employed to define the amount of the (de)refinement.

Expected result is drafted in figure below.

[[File:refinement_workflow.png|800px]]

== Implementation ==
The implementation of $hp$-adaptive solution procedure can be found [https://gitlab.com/e62Lab/public/2022_p_hp-adaptivity here].

== Solution procedure ==
=== Solve ===

First a numerical solution $\widehat{u}$ to the governing problem must be obtained. In general, the numerical treatment of a system of PDEs is done in several steps. First, the domain is discretized by positioning the nodes, then the linear differential operators in each computational node are approximated, and, finally, the system of PDEs is discretized and assembled into a large sparse linear system. To obtain a numerical solution $\widehat{u}$, the sparse system is solved.

=== Estimate (IMEX error indicator) ===
In the estimation step, critical areas with high error of the numerical solution are identified. Identifying such areas is not a trivial task. In rare cases, where a closed form solution to the governing system of PDEs exists, we can not only indicate such areas but also estimate the accuracy of the numerical solution. However, for real-world problems, which are often subject of numerical simulations, closed form solutions do not exist. Therefore, other objective metrics are needed to indicate the computational nodes with high error of the numerical solution. Numerical tools used in such cases are commonly referred to as ''error indicators''.

We have chosen to use our own error indication algorithm IMEX <ref name="imex"> Jančič, Mitja, Filip Strniša, and Gregor Kosec. "Implicit-Explicit Error Indicator based on Approximation Order." arXiv preprint arXiv:2204.01022 (2022).</ref> for implementation-related convenience reasons. This error indicator makes use of the implicitly obtained numerical solution and explicit operators (approximated by higher-order basis) to reconstruct the right-hand side of the governing problem.

To further explain basic idea of IMEX, let us define a PDE of type:
$$
\mathcal L u = f_{RHS},
$$
where $\mathcal L$ is a differential operator applied to the scalar field $u$ and $f_{RHS}$ is a function. To obtain an error indicator field $\eta$, the problem is first solved implicitly by using a lower order approximation of operators $\mathcal L$, $\mathcal L^{(lo)}_{im}$, obtaining the solution $u^{(im)}$ in the process. The implicitly computed field $u^{(im)}$ is then used to explicitly reconstruct the right-hand side of the governing problem, but this time using a higher order approximation of $\mathcal L$, $\mathcal L^{(hi)}_{ex}$, obtaining $f_{RHS}^{(ex)}$ in the process.
Finally, $f_{RHS}^{(ex)}$ and $f_{RHS}$ are compared $\eta = \left | f_{RHS} - f_{RHS}^{(ex)} \right |$ to indicate the error of the numerical solution.

=== Mark ===

After the error indicator $\eta$ had been obtained for each computational point in $\Omega$, a marking strategy is employed. The main objective of this step is to flag the nodes with too high or too low error indicator values in order to achieve a uniformly distributed precision of the numerical solution and reduce computational costs - by avoiding fine local field descriptions and high order approximations where this is not required. Furthermore, the marking strategy not only decides whether or not (de)refinement should take place at a given computational node, but also defines the type of the refinement procedure in cases where there are multiple to choose from.

In each adaptivity iteration, the marking strategy starts by checking the error indicator values $\eta _i$ for all computational nodes in the domain. If $\eta _i$ is greater than $\alpha \eta_{max}$ for a free model parameter $\alpha \in (0, 1)$, the node is marked for refinement. If $\eta _i$ is less than $\beta \eta_{max}$ for a free model parameter $\beta \in (0,1) \land \beta \leq \alpha$, the node is marked for derefinement. Otherwise, the node is left unmarked meaning no refinement should take place. The marking strategy can be summarized with a single equation
$$
\begin{cases}
\eta _i > \alpha \eta_{max}, & \text{ refine} \\
\beta \eta_{max} \leq \eta _i \leq \alpha \eta_{max}, & \text{ do nothing} \\
\eta _i < \beta \eta_{max}, & \text{ derefine}
\end{cases}.
$$

In the context of mesh-based methods, it has already been observed, that this strategy is indeed simple to implement, however, far from optimal in terms of achieving good convergence rates. Our studies have shown that a slight modification to the originally proposed strategy can significantly improve the performance of an $hp$-adaptive solution procedure in the context of mesh-free methods, but at the cost of higher number of free parameters - making it more difficult to understand. Nevertheless, the strategy is modified by introducing parameters $\left \{ \alpha_h, \beta_h \right \}$ and $\left \{ \alpha_p, \beta_p \right \}$ for separate treatment of $h$- and $p$-refinements respectively. A schematic is shown in figure below.

[[File:Screenshot from 2022-09-20 08-25-34.png|800px]]

=== Refine ===

After obtaining the list of nodes marked for modification, the refinement module is initialized. In this module, the local field description and local approximation order for the unmarked nodes are left unchanged, while the rest are further processed to determine other refinement-type-specific details - such as the amount of the (de)refinement.

We use the following $h$-refinement rule
$$
\label{eq:refinement}
h_i^{new}(\b p) = \frac{h_i^{old}}{\frac{\eta_i - \alpha \eta _{max}}{\eta_{max} - \alpha\eta_{max}}\big(\lambda - 1\big) + 1}
$$
for the dimensionless parameter $\lambda \in [1, \infty)$ allowing us to control the aggressiveness of the refinement - the larger the value, the greater the change. This refinement rule also conveniently refines the areas with higher error indicator values more than the ones that are closer to the upper refinement threshold $\alpha_h \eta_{max}$. Similarly, a derefinement rule is proposed
$$
\label{eq:derefinement}
h_i^{new}(\b p) = \frac{h_i^{old}}{\frac{\beta \eta _{max} - \eta_i}{\beta\eta_{max} - \eta_{min}}\big(\frac{1}{\vartheta} - 1\big) + 1},
$$
where parameter $\vartheta \in [1, \infty)$ allows us to control the aggressiveness of derefinement.

The same refinement and derefinement strategies were applied to control the local approximation order, except this time the value is rounded to the nearest integer. Similarly and for the same reasons as we did with the marking strategy, we consider a separate treatment of $h$- and $p$-adaptive procedures, by introducing (de)refinement aggressiveness parameters $\left \{\lambda_h, \vartheta_h \right \}$ and $\left \{\lambda_p, \vartheta_p \right \}$ for $h$- and $p$-refinement types respectively.

[[File:refinement_rules.png|800px]]

=== Finalization step ===

Before the 4 modules can be iteratively repeated, the domain is rediscretized taking into account the newly computed local internodal distances $h_i^{new}(\b p)$ and the local approximation orders $m_i^{new}(\b p)$. However, both are only known in the computational nodes, while global functions $\widehat{h}^{new}(\b p)$ and $\widehat{m}^{new}(\b p)$ over our entire domain space $\Omega$ are required. These are obtained by employing the Sheppard's inverse distance weighting interpolation.

Figure below schematically demonstrates 3 examples of $hp$-refinements. For the sake of demonstration, the refinement parameters for $h$- and $p$-adaptivity are the same, i.e. $\left \{ \alpha, \beta, \lambda, \vartheta \right \} = \left \{ \alpha_h, \beta_h, \lambda_h, \vartheta_h \right \} = \left \{ \alpha_p, \beta_p, \lambda_p, \vartheta_p \right \}$. Additionally, the derefinement aggressiveness $\vartheta$ and the upper threshold $\beta$ are kept constant, so, effectively, only the upper bound of refinement $\alpha$ and refinement aggressiveness $\lambda$ are altered. We observe that the effect refinement parameters is somewhat intuitive. The larger the aggressiveness $\lambda$ the better the local field description and the larger the number of nodes with high approximation order. Similar effect is observed with manipulating the upper refinement threshold $\alpha$, except the effect comes at a smoother manner. Also observe that all refined states were able to increase the accuracy of the numerical solution from the initial state.

[[File:refinement_demonstration.png|800px]]

== Examples ==
=== Example peak problem ===

The proposed $hp$-adaptive solution procedure is demonstrated on a synthetical example. We chose a 2-dimensional Poisson problem with exponentially strong source positioned at $\boldsymbol x_s = \Big (\frac{1}{2}, \frac{1}{3}\Big )$. This example is categorized as a difficult problem and is commonly used to test the performance of adaptive solution procedures. The problem has a tractable solution $u(\boldsymbol x)=e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$, which allows us to evaluate the precision of the numerical solution $\widehat{u}$, e.g.\ in terms of the infinity norm.

Governing equations are
$$
\nabla^2 u (\boldsymbol x) = 2a e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}(2a\left \| \boldsymbol x - \boldsymbol x_s \right \| - d) \quad \quad \text{in } \Omega, \\
u (\boldsymbol x) = e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2} \quad \quad \text{on } \Gamma_d, \\
\nabla u (\boldsymbol x) = -2a(\boldsymbol x - \boldsymbol x_s)e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2} \quad \quad\text{on } \Gamma_n,
$$

Figure below shows example indicator fields for initial iteration, intermediate iteration and iteration that achieved the best accuracy of the numerical solution. The third column shows the IMEX error indicator. We can see that the IMEX successfully located the position of the strong source at $\boldsymbol x_s = \Big (\frac{1}{2}, \frac{1}{3}\Big )$ as the highest indicator values are seen in its neighborhood. Moreover, the second column shows that the accuracy of the numerical solution and uniformity of error distribution have both been significantly improved by the $hp$-adaptive solution procedure, further proving that the IMEX can be successfully employed as a reliable error indicator.

[[File:refinement_demonstration_2d.png|800px]]

The behavior of IMEX over 70 adaptivity iterations is further studied in figure below. We are pleased to see that the convergence limit of the indicator around iteration $N_{iter}=60$ agrees well with the convergence limit of the numerical solution.

[[File:error_indicator_convergence.png|800px]]

Finally, the convergence behavior of the proposed $hp$-adaptive solution procedure is studied. In addition to the convergence of a single $hp$-adaptive run, the convergences obtained without the use of refinement procedures, i.e. solutions obtained with uniform internodal spacing and approximation orders over the entire domain, are shown in figure below. The figure evidently shows that a $hp$-adaptive solution procedure was able to notably improve the numerical solution in terms of accuracy and required computational points.

[[File:convergence_of_h_hp.png|800px]]

=== Application to linear elasticity problems: Fretting fatigue contact ===

The application of the proposed $hp$-adaptive solution procedure is further expanded to study a linear elasticity problem. Specifically, we obtain a $hp$-refined solution to fretting fatigue contact problem <ref name="fwo">Pereira, Kyvia, et al. "On the convergence of stresses in fretting fatigue." Materials 9.8 (2016): 639.</ref>, for which the authors believe no closed form solution is known.

The problem formulation is given here <ref name="slak">Slak, Jure, and Gregor Kosec. "Adaptive radial basis function–generated finite differences method for contact problems." International Journal for Numerical Methods in Engineering 119.7 (2019): 661-686.</ref>.

Figure below shows an example of $hp$-refined solution to fretting fatigue problem at the last adaptivity iteration with $N=46\,626$ computational nodes. We see that the solution procedure successfully located the two critical points, i.e. the fixed top left corner with a stress singularity and the area in the middle of the top edge where the contact is simulated. Note that the highest stress values (approximately $2$ times higher) have been computed in the singularity in the top left corner, but those nodes are not shown since our focus is shifted towards the area under the contact.

[[File:fwo_solution_hp.png|800px]]

For a detailed analysis, we look at the surface traction $\sigma_{xx}$, as it is often used to determine the location of crack initiation. The surface traction is shown in figure below for $6$ selected adaptivity iterations. The mesh-free nodes are colored according to the local approximation order enforced by the $hp$-adaptive solution procedure. The message of this figure is twofold. Firstly, it is clear that the proposed IMEX error indicator can be successfully used in linear elasticity problems, and secondly, we observe that the $hp$-adaptive solution procedure successfully approximated the surface traction in the neighborhood of the contact. In the process, the local field description under the contact has been significantly improved and also the local approximation orders have taken a non-trivial distribution.

[[File:fwo_abaqus.png|800px]]

The surface traction is additionally accompanied with the FEM results on a much denser mesh with more than 100\,000 DOFs obtained with the commercial solver Abaqus\textsuperscript{\textregistered}. To compute the absolute difference between the two methods, the mesh-free solution has been interpolated to Abaqus's computational points using the Sheppard's inverse distance weighting interpolation with $2$ closest neighbors. We see that the absolute difference is decreasing with the number of adaptivity iterations, finally settling at approximately 2 % of the maximum value under the contact at initial iteration. The highest absolute difference is expectedly located at the edges of the contact, i.e. around $x=a$ and $x=-a$, while the difference in the rest of the contact is even smaller.

=References=
<references/>

Go back to [[Medusa#Examples|Examples]].

HP-adaptivity

2022-09-20T07:47:57Z

Mitjajancic:

Adaptive solution procedures are essential in problems where the accuracy of the numerical solution varies spatially and are currently subject of intensive studies. Two conceptually different adaptive approaches have been proposed, namely $p$-adaptivity or $h$-, $r$-adaptivity. In $p$-adaptivity, the accuracy of the numerical solution is varied by changing the order of approximation, while in $h$- and $r$-adaptivity the resolution of the spatial discretization is adjusted for the same purpose. In the $h$-adaptive approach, nodes are added or removed from the domain as needed, while in the $r$-adaptive approach, the total number of nodes remains constant - the nodes are only repositioned with respect to the desired accuracy. Ultimately, $h$- and $p$-adaptivities can be combined to form the so-called $hp$-adaptivity, where the accuracy of the solution is controlled with the order of the method and the resolution of the spatial discretization.

__FORCETOC__

== Basic concept ==

The proposed $hp$-adaptive solution procedure follows the well-established paradigm based on an iterative loop, where each iteration step consists of four modules:
* '''Solve''' - A numerical solution $\widehat{u}$ is obtained.
* '''Estimate''' - An estimate of the spatial accuracy of the numerical solution is calculated using error indicators.
* '''Mark''' - Depending on the error indicator values $\eta _i$, a marking strategy is used to mark the computational nodes for (de)refinement.
* '''Refine''' - Refinement strategy is employed to define the amount of the (de)refinement.

Expected result is drafted in figure below.

[[File:refinement_workflow.png|800px]]

== Implementation ==
The implementation of $hp$-adaptive solution procedure can be found [https://gitlab.com/e62Lab/public/2022_p_hp-adaptivity here].

== Solution procedure ==
=== Solve ===

First a numerical solution $\widehat{u}$ to the governing problem must be obtained. In general, the numerical treatment of a system of PDEs is done in several steps. First, the domain is discretized by positioning the nodes, then the linear differential operators in each computational node are approximated, and, finally, the system of PDEs is discretized and assembled into a large sparse linear system. To obtain a numerical solution $\widehat{u}$, the sparse system is solved.

=== Estimate (IMEX error indicator) ===
In the estimation step, critical areas with high error of the numerical solution are identified. Identifying such areas is not a trivial task. In rare cases, where a closed form solution to the governing system of PDEs exists, we can not only indicate such areas but also estimate the accuracy of the numerical solution. However, for real-world problems, which are often subject of numerical simulations, closed form solutions do not exist. Therefore, other objective metrics are needed to indicate the computational nodes with high error of the numerical solution. Numerical tools used in such cases are commonly referred to as ''error indicators''.

We have chosen to use our own error indication algorithm IMEX <ref name="imex"> Jančič, Mitja, Filip Strniša, and Gregor Kosec. "Implicit-Explicit Error Indicator based on Approximation Order." arXiv preprint arXiv:2204.01022 (2022).</ref> for implementation-related convenience reasons. This error indicator makes use of the implicitly obtained numerical solution and explicit operators (approximated by higher-order basis) to reconstruct the right-hand side of the governing problem.

To further explain basic idea of IMEX, let us define a PDE of type:
$$
\mathcal L u = f_{RHS},
$$
where $\mathcal L$ is a differential operator applied to the scalar field $u$ and $f_{RHS}$ is a function. To obtain an error indicator field $\eta$, the problem is first solved implicitly by using a lower order approximation of operators $\mathcal L$, $\mathcal L^{(lo)}_{im}$, obtaining the solution $u^{(im)}$ in the process. The implicitly computed field $u^{(im)}$ is then used to explicitly reconstruct the right-hand side of the problem~\eqref{eq: general PDE}, but this time using a higher order approximation of $\mathcal L$, $\mathcal L^{(hi)}_{ex}$, obtaining $f_{RHS}^{(ex)}$ in the process.
Finally, $f_{RHS}^{(ex)}$ and $f_{RHS}$ are compared $\eta = \left | f_{RHS} - f_{RHS}^{(ex)} \right |$ to indicate the error of the numerical solution.

=== Mark ===

After the error indicator $\eta$ had been obtained for each computational point in $\Omega$, a marking strategy is employed. The main objective of this step is to flag the nodes with too high or too low error indicator values in order to achieve a uniformly distributed precision of the numerical solution and reduce computational costs - by avoiding fine local field descriptions and high order approximations where this is not required. Furthermore, the marking strategy not only decides whether or not (de)refinement should take place at a given computational node, but also defines the type of the refinement procedure in cases where there are multiple to choose from.

In each adaptivity iteration, the marking strategy starts by checking the error indicator values $\eta _i$ for all computational nodes in the domain. If $\eta _i$ is greater than $\alpha \eta_{max}$ for a free model parameter $\alpha \in (0, 1)$, the node is marked for refinement. If $\eta _i$ is less than $\beta \eta_{max}$ for a free model parameter $\beta \in (0,1) \land \beta \leq \alpha$, the node is marked for derefinement. Otherwise, the node is left unmarked meaning no refinement should take place. The marking strategy can be summarized with a single equation
$$
\begin{cases}
\eta _i > \alpha \eta_{max}, & \text{ refine} \\
\beta \eta_{max} \leq \eta _i \leq \alpha \eta_{max}, & \text{ do nothing} \\
\eta _i < \beta \eta_{max}, & \text{ derefine}
\end{cases}.
$$

In the context of mesh-based methods, it has already been observed, that this strategy is indeed simple to implement, however, far from optimal in terms of achieving good convergence rates. Our studies have shown that a slight modification to the originally proposed strategy can significantly improve the performance of an $hp$-adaptive solution procedure in the context of mesh-free methods, but at the cost of higher number of free parameters - making it more difficult to understand. Nevertheless, the strategy is modified by introducing parameters $\left \{ \alpha_h, \beta_h \right \}$ and $\left \{ \alpha_p, \beta_p \right \}$ for separate treatment of $h$- and $p$-refinements respectively. A schematic is shown in figure below.

[[File:Screenshot from 2022-09-20 08-25-34.png|800px]]

=== Refine ===

After obtaining the list of nodes marked for modification, the refinement module is initialized. In this module, the local field description and local approximation order for the unmarked nodes are left unchanged, while the rest are further processed to determine other refinement-type-specific details - such as the amount of the (de)refinement.

We use the following $h$-refinement rule
$$
\label{eq:refinement}
h_i^{new}(\b p) = \frac{h_i^{old}}{\frac{\eta_i - \alpha \eta _{max}}{\eta_{max} - \alpha\eta_{max}}\big(\lambda - 1\big) + 1}
$$
for the dimensionless parameter $\lambda \in [1, \infty)$ allowing us to control the aggressiveness of the refinement - the larger the value, the greater the change. This refinement rule also conveniently refines the areas with higher error indicator values more than the ones that are closer to the upper refinement threshold $\alpha_h \eta_{max}$. Similarly, a derefinement rule is proposed
$$
\label{eq:derefinement}
h_i^{new}(\b p) = \frac{h_i^{old}}{\frac{\beta \eta _{max} - \eta_i}{\beta\eta_{max} - \eta_{min}}\big(\frac{1}{\vartheta} - 1\big) + 1},
$$
where parameter $\vartheta \in [1, \infty)$ allows us to control the aggressiveness of derefinement.

The same refinement and derefinement strategies were applied to control the local approximation order, except this time the value is rounded to the nearest integer. Similarly and for the same reasons as we did with the marking strategy, we consider a separate treatment of $h$- and $p$-adaptive procedures, by introducing (de)refinement aggressiveness parameters $\left \{\lambda_h, \vartheta_h \right \}$ and $\left \{\lambda_p, \vartheta_p \right \}$ for $h$- and $p$-refinement types respectively.

[[File:refinement_rules.png|800px]]

=== Finalization step ===

Before the 4 modules can be iteratively repeated, the domain is rediscretized taking into account the newly computed local internodal distances $h_i^{new}(\b p)$ and the local approximation orders $m_i^{new}(\b p)$. However, both are only known in the computational nodes, while global functions $\widehat{h}^{new}(\b p)$ and $\widehat{m}^{new}(\b p)$ over our entire domain space $\Omega$ are required. These are obtained by employing the Sheppard's inverse distance weighting interpolation.

Figure below schematically demonstrates 3 examples of $hp$-refinements. For the sake of demonstration, the refinement parameters for $h$- and $p$-adaptivity are the same, i.e.\ $\left \{ \alpha, \beta, \lambda, \vartheta \right \} = \left \{ \alpha_h, \beta_h, \lambda_h, \vartheta_h \right \} = \left \{ \alpha_p, \beta_p, \lambda_p, \vartheta_p \right \}$. Additionally, the derefinement aggressiveness $\vartheta$ and the upper threshold $\beta$ are kept constant, so, effectively, only the upper bound of refinement $\alpha$ and refinement aggressiveness $\lambda$ are altered. We observe that the effect refinement parameters is somewhat intuitive. The larger the aggressiveness $\lambda$ the better the local field description and the larger the number of nodes with high approximation order. Similar effect is observed with manipulating the upper refinement threshold $\alpha$, except the effect comes at a smoother manner. Also observe that all refined states were able to increase the accuracy of the numerical solution from the initial state.

[[File:refinement_demonstration.png|800px]]

== Examples ==
=== Example peak problem ===

The proposed $hp$-adaptive solution procedure is demonstrated on a synthetical example. We chose a 2-dimensional Poisson problem with exponentially strong source positioned at $\boldsymbol x_s = \Big (\frac{1}{2}, \frac{1}{3}\Big )$. This example is categorized as a difficult problem and is commonly used to test the performance of adaptive solution procedures. The problem has a tractable solution $u(\boldsymbol x)=e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}$, which allows us to evaluate the precision of the numerical solution $\widehat{u}$, e.g.\ in terms of the infinity norm.

Governing equations are
$$
\nabla^2 u (\boldsymbol x) = 2a e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2}(2a\left \| \boldsymbol x - \boldsymbol x_s \right \| - d) \quad \quad \text{in } \Omega, \\
u (\boldsymbol x) = e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2} \quad \quad \text{on } \Gamma_d, \\
\nabla u (\boldsymbol x) = -2a(\boldsymbol x - \boldsymbol x_s)e^{-a \left \| \boldsymbol x - \boldsymbol x_s \right \|^2} \quad \quad\text{on } \Gamma_n,
$$

Figure below shows example indicator fields for initial iteration, intermediate iteration and iteration that achieved the best accuracy of the numerical solution. The third column shows the IMEX error indicator. We can see that the IMEX successfully located the position of the strong source at $\boldsymbol x_s = \Big (\frac{1}{2}, \frac{1}{3}\Big )$ as the highest indicator values are seen in its neighborhood. Moreover, the second column shows that the accuracy of the numerical solution and uniformity of error distribution have both been significantly improved by the $hp$-adaptive solution procedure, further proving that the IMEX can be successfully employed as a reliable error indicator.

[[File:refinement_demonstration_2d.png|800px]]

The behavior of IMEX over 70 adaptivity iterations is further studied in figure below. We are pleased to see that the convergence limit of the indicator around iteration $N_{iter}=60$ agrees well with the convergence limit of the numerical solution.

[[File:error_indicator_convergence.png|800px]]

Finally, the convergence behavior of the proposed $hp$-adaptive solution procedure is studied. In addition to the convergence of a single $hp$-adaptive run, the convergences obtained without the use of refinement procedures, i.e. solutions obtained with uniform internodal spacing and approximation orders over the entire domain, are shown in figure below. The figure evidently shows that a $hp$-adaptive solution procedure was able to notably improve the numerical solution in terms of accuracy and required computational points.

[[File:convergence_of_h_hp.png|800px]]

=== Application to linear elasticity problems: Fretting fatigue contact ===

The application of the proposed $hp$-adaptive solution procedure is further expanded to study a linear elasticity problem. Specifically, we obtain a $hp$-refined solution to fretting fatigue contact problem <ref name="fwo">Pereira, Kyvia, et al. "On the convergence of stresses in fretting fatigue." Materials 9.8 (2016): 639.</ref>, for which the authors believe no closed form solution is known.

The problem formulation is given here <ref name="slak">Slak, Jure, and Gregor Kosec. "Adaptive radial basis function–generated finite differences method for contact problems." International Journal for Numerical Methods in Engineering 119.7 (2019): 661-686.</ref>.

Figure below shows an example of $hp$-refined solution to fretting fatigue problem at the last adaptivity iteration with $N=46\,626$ computational nodes. We see that the solution procedure successfully located the two critical points, i.e. the fixed top left corner with a stress singularity and the area in the middle of the top edge where the contact is simulated. Note that the highest stress values (approximately $2$ times higher) have been computed in the singularity in the top left corner, but those nodes are not shown since our focus is shifted towards the area under the contact.

[[File:fwo_solution_hp.png|800px]]

For a detailed analysis, we look at the surface traction $\sigma_{xx}$, as it is often used to determine the location of crack initiation. The surface traction is shown in Figure~\ref{fig:fwo_conv} for $6$ selected adaptivity iterations. The mesh-free nodes are colored according to the local approximation order enforced by the $hp$-adaptive solution procedure. The message of this figure is twofold. Firstly, it is clear that the proposed IMEX error indicator can be successfully used in linear elasticity problems, and secondly, we observe that the $hp$-adaptive solution procedure successfully approximated the surface traction in the neighborhood of the contact. In the process, the local field description under the contact has been significantly improved and also the local approximation orders have taken a non-trivial distribution.

[[File:fwo_abaqus.png|800px]]

The surface traction in Figure~\ref{fig:fwo_conv} is additionally accompanied with the FEM results on a much denser mesh with more than 100\,000 DOFs obtained with the commercial solver Abaqus\textsuperscript{\textregistered}. To compute the absolute difference between the two methods, the mesh-free solution has been interpolated to Abaqus's computational points using the Sheppard's inverse distance weighting interpolation with $2$ closest neighbors. We see that the absolute difference is decreasing with the number of adaptivity iterations, finally settling at approximately 2 \% of the maximum value under the contact at initial iteration. The highest absolute difference is expectedly located at the edges of the contact, i.e.\ around $x=a$ and $x=-a$, while the difference in the rest of the contact is even smaller.

=References=
<references/>

Go back to [[Medusa#Examples|Examples]].

HP-adaptivity

2022-09-20T07:06:11Z

Mitjajancic:

TODO

File:Fwo abaqus.png

2022-09-20T06:31:12Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Error indicator convergence.png

2022-09-20T06:31:11Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Convergence of h hp.png

2022-09-20T06:31:11Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Fwo solution hp.png

2022-09-20T06:31:11Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Refinement demonstration.png

2022-09-20T06:31:10Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Refinement demonstration 2d.png

2022-09-20T06:31:10Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Screenshot from 2022-09-20 08-25-34.png

2022-09-20T06:31:09Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Refinement rules.png

2022-09-20T06:31:09Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

File:Refinement workflow.png

2022-09-20T06:24:22Z

Mitjajancic: File uploaded with MsUpload

File uploaded with MsUpload

Adaptivity

2022-09-20T05:48:34Z

Mitjajancic: typo

Go back to [[Medusa#Examples|Examples]].

Solutions to many physical problems governed by partial differential equations (PDE) often significantly vary in magnitude throughout the problem domain. Although in some special cases the areas with high error are known in advance, in general the error distribution is unknown beforehand. Adaptive techniques for solving PDEs are a standard way of dealing with this problem, where problematic regions are iteratively refined. A step further is automatic adaptivity, where problematic regions are chosen automatically using an error indicator and then refined, until certain error threshold is reached. Below, we show some examples of fully automatic adaptivity in Medusa.

__FORCETOC__

== Basic concept ==

The adaptive methodology in this paper behaves similarly to "remeshing" used commonly in FEM.
Some initial (possibly variable) nodal spacing $h^0$ is chosen, as well as its lower and upper bounds $h_L$ and $h_U$, respectively.
Domain $\Omega$ is filled with nodes, conforming to $h^0$ and the solution $u^0$ is obtained.
An error indicator is employed to determine which nodes should be (de)refined and the nodal density $h^0$ is altered appropriately.
This adaptive cycle below is repeated until the convergence criterion is met. The procedure on $j$-th iteration is written in more detail below:

<ol>
<li>Fill $\Omega$ with nodes conforming to $h^j$.</li>
<li>Solve the problem to obtain $u^j$.</li>
<li>Compute the error indicator values $\varepsilon_i^j$ for each node $p_i$.</li>
<li>If the mean of $\varepsilon_i^j$ is below some tolerance $\varepsilon$ return $u^j$ as the solution and stop.</li>
<li>Adapt $h^j$ to obtain $h^{j+1}$.</li>
</ol>

More details can be found in our paper <ref name="SlakAdaptive">Slak, J., & Kosec, G. Adaptive radial basis function‐generated finite differences method for contact problems. International Journal for Numerical Methods in Engineering, 2019.</ref>.

== Node density adaptation ==

The existing nodal spacing function $h^j$ is evaluated at nodes $p_i$ to obtain values $h_{i,j} = h^j(p_i)$ . These values $h_{i,j}$ are modified by a density factor $f_i$ as
$$
h_{i,j+1} = \min(\max(h_{i,j} / f_i, h_L(p_i)), h_U(p_i)),
$$
where density factor $f_i$ is computed as
$$
f_i = \begin{cases}
1 + \frac{\eta - \varepsilon_i}{\eta - m} (\frac{1}{\beta} - 1), & \varepsilon_i \leq \eta, \quad \text{i.e. decrease the density} \\
1, & \eta < \varepsilon_i < \varepsilon, \quad \text{i.e. no change in density}\\
1 + \frac{\varepsilon_i - \varepsilon}{M - \varepsilon} (\alpha - 1), & \varepsilon_i \geq \varepsilon, \quad \text{i.e. increase the density}
\end{cases}
$$
and $\alpha$ represents the refine aggressiveness, $\beta$ the derefine aggressiveness, $\varepsilon$ the refinement threshold, $\eta$ the derenfinement threshold,
and $m = \min_i \varepsilon_i$ is the minimal and $M = \max_i \varepsilon_i$ is the maximal value of the error indicator.
Note that setting $\alpha=1$ or $\beta=1$ disables refinement and derefinement, respectively.

This adaptation is illustrated in the figure below.

[[File:density_change.png|1193px]]

== Error indicators ==
The work on error indicators is ongoing. For now, we use an ad hoc error indicator
$$\varepsilon_i = \operatorname{std}_{j \in I_i}(u_j),$$
which represents the standard deviation of function values over all stencil nodes of a given node $p_i$.

== Numerical exmaples ==
Below are several numerical examples where adaptivity has been tested or used to obtain solutions.

The errors $e_1$, $e_2$, $e_\infty$ and $e_E$ refer to relative discrete $L^1$, $L^2$, $L_\infty$ and energy norm errors, respectively. These are evaluated in the computation nodes or on a denser grid by reinterpolation.

Relative node density $\rho_i$ is calculated by computing the distances $d_i$ from each node to its closest neighbour and expressing them as the logarithmic ratio against the maximal distance: $\rho_i = -\log_2(d_i / \max_j d_j)$.

=== L shaped domain ===

The L shaped domain problem is defined on $\Omega = [-1, 1]^2 \setminus [0, 1] \times [-1, 0]$. The Laplace problem $\nabla^2u = 0$ with the solution $u = r^{\frac{2}{3}} \sin(\frac{2}{3}\theta)$ given in polar coordinates.

BF-FD method with Polyharmonic splines augmented with monomials up to and including 2nd order was used to approximate the differential operators. The stencils for each node were chosen by simply selecting the closest $n=15$ nodes. The resulting sparse system was solved using the Intel ® MKL Pardiso sparse solver. Both uniform and fully adaptive refinement was tested. The adaptive procedure was run with $\alpha=3$, $\varepsilon = 10^{-2}$, $\beta=1$ and $\eta=0$.

The errors under uniform (left) and adaptive (right) refinement are shown below.

[[File:L_shape_uniform_error.png|500px]][[File:L_shape_adaptive_error.png|500px]]

The error (left) and the nodal density (right) during the adaptive iteration are shown below.

[[File:L_shape_progress.png|600px]]

=== Disk under stress ===

The next case is disk under pressure case from [[Linear elasticity]] solving the Caucy-Navier equation of [[Solid Mechanics]].
The problem considers one quarter of the disk illustrated below by the domain $\Omega$. The spacing $\gamma$ represents the distance from the singularity.

[[File:disk.png|350px]]

Computed $\sigma_{xy}$ stress profiles for three considered cases. The lower three images show close-ups of the computed peaks.

[[File:disk_profiles.png|1050px]]

Error of the numerical solution of the compressed disk problem during adaptive iteration.

[[File:disk_error.png|977px]]

Error indicator value and node density in the last iteration.

[[File:disk_err_den.png|695px]]

Comparison of uniform and adaptive refinement.

[[File:disk_uni_ada_err.png|971px]]

=== Hertzian contact ===

See [[Hertzian_contact]] for the description.

Errors and node counts during adaptive iteration for the solution of the Hertzian problem. Observe the initial decrease due to derefinement.

[[File:her_err_n.png|921px]]

Top surface stress profiles during the iteration. Only even iterations are shown for brevity.

[[File:her_profile.png|840px]]

Comparing adaptive (left) and manual (right) nodal distributions densities for the Hertzian problem. The manual density was used in the refinement paper<ref name="SlakRef">Slak, J., & Kosec, G. (2019). Refined Meshless Local Strong Form solution of Cauchy–Navier equation on an irregular domain. Engineering analysis with boundary elements, 100, 3-13.</ref>.

[[File:her_density.png|906px]]

=== Fretting fatigue contact ===

See [[Fretting fatigue case]] for the description.

Surface traction $\sigma_{xx}$ under contact in four subsequent adaptive iterations during the solution of the fretting fatigue problem.

[[File:fwo_profiles.png|627px]]

Nodes of four subsequent adaptive iterations during the solution of the fretting fatigue problem, coloured by von Mises stress.

[[File:fwo_iters.png|819px]]

Error of surface traction using two reference solutions.

[[File:fwo_error.png|505px]]

=== Bousinesq's problem ===

See [[Linear_elasticity#Point_contact_3D]] for more details.

Error with respect to the number of nodes in the adaptive iteration of the solution of the Boussinesq's problem.
The right axis shows the number of computational nodes.

[[File:bou_err.png|778px]]

Stress profiles along the body diagonal from $[-1, -1, -1]$ to $[-\gamma, -\gamma, -\gamma]$ in adaptive iteration when solving the
Boussinesq's problem.

[[File:bou_profile.png|1063px]]

The obtained solution in the final iteration with an enlarged portion around the contact area. Both solutions are
plotted only in the computation nodes to also show the final nodal distribution. Nodes are coloured proportional to computed
values of von Mises stress.

[[File:bou_sol.png|1165px]]

=== Helmholz equation ===

The equation $\nabla^2 u + \frac{1}{(\alpha+r)^4} u = f$ with the solution $u = sin(\frac{1}{\alpha+r})$ is solved adaptively. The value $\alpha = \frac{1}{5\pi}$ was used and the right hand side $f$ was computed from $u$.

Error under uniform (left) and adaptive (right) refinement.

[[File:osc2d_uniform_error.png|500px]][[File:osc2d_adaptive_error.png|500px]]

Solution and the nodal density in the final iteration.

[[File:osc2d_adaptive_sol.png|500px]][[File:osc2d_adaptive_den.png|500px]]

Error under adaptive refinement for the 3D equation.

[[File:osc3d_adaptive_error.png|500px]]

Solution and the nodal density in the final iteration of the 3D solution.

[[File:osc3d_adaptive_sol.png|500px]][[File:osc3d_adaptive_den.png|500px]]

=References=
<references/>

Go back to [[Medusa#Examples|Examples]].

Medusa

2022-09-20T05:41:55Z

Mitjajancic: Added hp-adaptivity page.

'''Welcome to the Medusa wiki. To visit the main website, go to [http://e6.ijs.si/medusa/ http://e6.ijs.si/medusa/].'''

In [http://e6.ijs.si/ParallelAndDistributedSystems/ Parallel and Distributed Systems Laboratory] we are working on a C++ library that is first and foremost focused on tools for solving Partial Differential Equations by meshless methods. The basic idea is to create generic codes for tools that are needed for solving not only PDEs but many other problems, e.g. Moving Least Squares approximation, $k$-d tree, domain generation engines, etc.
We call this open source meshless project [http://e6.ijs.si/medusa/ Medusa: Coordinate Free Meshless Method implementation (MM)].

Technical details about code and examples can be found on our [http://e6.ijs.si/medusa/docs/ documentation page] and [https://gitlab.com/e62Lab/medusa Gitlab repository]. [[File:C.png|100px||link=https://gitlab.com/e62Lab/medusa|alt=Alt text|code]] [[File:doxygen.png|100px|link=http://e6.ijs.si/medusa/docs/|alt=Alt text|Documentation page]]

This wiki site is meant for more relaxed discussions about general principles, possible and already implemented applications, preliminary analyses, etc.
Note, that there are many grammatical mistakes, typos, stupid sentences, etc. This wiki is meant for quick information exchange and therefore we do not invest a lot of energy into styling :).

== Documentation ==
* [https://gitlab.com/e62Lab/medusa Code on Gitlab]
* [[How to build | Installation and building]]
* [[Including this library in your project | Including this library in your project]]
* [[Testing | Running tests]]
* [http://e6.ijs.si/medusa/docs/ Technical documentation]
* [[Coding style | Coding style]]
* [[Wiki editing guide | Wiki editing and backup guide]]

== Examples ==
In this section we present exact examples. Each of the below solutions can be found also in in the repository under examples. More explanation about the physical background and solution procedure can be found in following sections.
* [[Philosophy of examples and how to run them]]
* [[Poisson's equation]]
* [[Heat equation]]
* [[Linear elasticity]]
* [[Complex-valued problems]]
* [[Coupled domains]]
* [[Parametric domains | Parametric domains – Curved surface with variable density]]
* [[NURBS domains | Domains modeled with non-uniform rational basis spline's (NURBS)]]
* [[Determining the interior of the domain by oversampling the boundary]]
* [[Computer-aided design - Importing IGES and STEP files]]
* [[Realistic 3D models]]
* [[Customization]]
* [[Ghost nodes]]
* [[Electromagnetic scattering]]
* [[Schrödinger equation]]
* [[Wave equation]]
* [[Cahn-Hilliard equation]]
* [[Non-Newtonian_fluid_example | Non-Newtonian fluid]]
* [[Meshless Lattice Boltzmann method]]

== Building blocks ==
Medusa is modular coordinate-free parallel implementation of a numerical framework designed, but not limited to, for solving PDEs. In this section we present main modules of the library that can be also used as a standalone tools.
* [[Positioning of computational nodes]]
* [[Relaxation of the nodal distribution]]
* [[Refinement of the nodal distribution]]
* [[k-d tree|''k''-d tree]] and other spatial search structures
* Solving [[Solving linear systems | linear systems]], [[Solving overdetermined systems | overdetermined]] and [[Solving underdetermined systems | underdetermined]]
* [[Weighted Least Squares (WLS)]]
* [[Computation of shape functions]]
* [[Meshless Local Strong Form Method (MLSM)]]
* [[Radial basis function-generated finite differences (RBF-FD)]]
* [[Ghost nodes (theory)]]
* [[Integrators for time stepping]]
* [[RBF Interpolation]]

== Discussions / Applications ==
This section is meant for general discussion about the physical background of the examples, the solution procedures, various applications, etc. Note, that code snippets presented in discussion might not reflect the actual state of Medusa.
* [[Convection Diffusion equation | Convection Diffusion equation]]
* [[Adaptivity]]
** [[HP-adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
** [[Plasticity]]
* [[Fluid Mechanics]]
** [[Lid driven cavity]]
** [[de Vahl Davis natural convection test]]
** [[Natural convection in 3D irregular domain]]
** [[Natural convection from heated cylinder]]
** [[Natural convection between concentric cylinders]]
** [[Non-Newtonian fluid]]
* [[Computational electromagnetics]]
** [[Triple dielectric step in 1D]]
** [[Scattering from an infinite cylinder]]
** [[Point source near an anisotropic lens]]
* Other applications
** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]
** [[Heart rate variability detection]]
** [[Bioheat equation]]

== Performance analyses ==
* [[Execution on Intel® Xeon Phi™ co-processor]]
* [[1D MLSM and FDM comparison]]
* [[:File:tech_report.pdf|Execution overheads due to clumsy types::technical report]] [[File:pdf-file.gif]]
* [[Solving sparse systems]]
* [[Eigen paralelization]]

== Last changes ==
<news unique=1 limit = 5>
*{{{timeanddate}}} :: {{{title}}}

</news>

== Miscellaneous ==
* FAQ - [[Frequently asked questions]].
* [[List of wiki contributors]]
* List of library contributors: [http://e6.ijs.si/medusa/about#about-contributors See the official website]

== References ==

For pre-prints check https://e6.ijs.si/ParallelAndDistributedSystems/publications/

* J. Slak, G. Kosec; Medusa : A C++ library for solving PDEs using strong form mesh-free methods, ACM transactions on mathematical software, vol. 47, 2021 [DOI: 10.1145/3450966]
* M. Depolli, J. Slak, G. Kosec; Parallel domain discretization algorithm for RBF-FD and other meshless numerical methods for solving PDEs, Computers & Structures, 2022 [DOI: 10.1016/j.compstruc.2022.106773]
* U. Duh, G. Kosec, J. Slak; Fast variable density node generation on parametric surfaces with application to mesh-free methods, SIAM journal on scientific computing, vol. 43, 2021 [DOI: 10.1137/20M1325642]
* M. Jančič, J. Slak, G. Kosec; Monomial augmentation guidelines for RBF-FD from accuracy versus computational time perspective, Journal of scientific computing, vol. 87, 2021 [DOI: 10.1007/s10915-020-01401-y]
* Slak J., Kosec G., On Generation of Node Distributions for Meshless PDE Discretizations, SIAM Journal on Scientific Computing, 41(5), A3202
* Slak J., Kosec G. Adaptive radial basis function-generated finite differences method for contact problems. International journal for numerical methods in engineering, ISSN 0029-5981 [http://www-e6.ijs.si/ParallelAndDistributedSystems/pdf/32230439.pdf manuscript]
* Slak J., Kosec G.; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain. Engineering analysis with boundary elements. 2018;11 ; [http://comms.ijs.si/~gkosec/data/papers/31107623.pdf manuscript]
* Depolli, M., Kosec, G., Assessment of differential evolution for multi-objective optimization in a natural convection problem solved by a local meshless method. Engineering optimization, 2017, vol. 49, no. 4, pp. 675-692 ;[http://comms.ijs.si/~gkosec/data/papers/29639719.pdf manuscript]
* Kosec G., A local numerical solution of a fluid-flow problem on an irregular domain. Advances in engineering software. 2016;7 ; [29512743] ; [http://comms.ijs.si/~gkosec/data/papers/29512743.pdf manuscript]
* Kosec G., Trobec R., Simulation of semiconductor devices with a local numerical approach. Engineering analysis with boundary elements. 2015;69-75; [27912487] ; [http://comms.ijs.si/~gkosec/data/papers/27912487.pdf manuscript]
* Kosec G., Šarler B., Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method. Engineering analysis with boundary elements. 2014;36-44; [http://comms.ijs.si/~gkosec/data/papers/3218939.pdf manuscript]
* Kosec G., Depolli M., Rashkovska A., Trobec R., Super linear speedup in a local parallel meshless solution of thermo-fluid problem. Computers & Structures. 2014;133:30-38; [http://comms.ijs.si/~gkosec/data/papers/27339815.pdf manuscript]
* Kosec G., Zinterhof P., Local strong form meshless method on multiple Graphics Processing Units. Computer modeling in engineering & sciences. 2013;91:377-396; [http://comms.ijs.si/~gkosec/data/papers/26785063.pdf manuscript]
* Kosec G., Šarler B., H-adaptive local radial basis function collocation meshless method. Computers, materials & continua. 2011;26:227-253; [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerBurgers.pdf manuscript]
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321; [http://comms.ijs.si/~gkosec/data/papers/EABE2499.pdf manuscript]
* Kosec G, Zaloznik M, Sarler B, Combeau H. A Meshless Approach Towards Solution of Macrosegregation Phenomena. CMC: Computers, Materials, & Continua. 2011;580:1-27 [http://comms.ijs.si/~gkosec/data/papers/KosecZaloznikSarlerCombeauSegregation.pdf manuscript]
* Kosec G, Sarler B. Solution of thermo-fluid problems by collocation with local pressure correction. International Journal of Numerical Methods for Heat & Fluid Flow. 2008;18:868-82 [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerNS2008.pdf manuscript]
* Trobec R., Kosec G., Parallel Scientific Computing, ISBN: 978-3-319-17072-5 (Print) 978-3-319-17073-2.
* Slak, J., Kosec, G.. Detection of heart rate variability from a wearable differential ECG device., MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938, pp 450-455.
* Kolman, M., Kosec, G. Correlation between attenuation of 20 GHz satellite communication link and liquid water content in the atmosphere. MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938. pp. 308-313.

==Related pages==
* http://e6.ijs.si/ParallelAndDistributedSystems/products/medusa/

HP-adaptivity

2022-09-20T05:41:13Z

Mitjajancic: Mitjajancic moved page Hp-adaptivity to HP-adaptivity

Test

Hp-adaptivity

2022-09-20T05:41:13Z

Mitjajancic: Mitjajancic moved page Hp-adaptivity to HP-adaptivity

#REDIRECT [[HP-adaptivity]]

HP-adaptivity

2022-09-20T05:40:15Z

Mitjajancic: Mitjajancic moved page $hp$-adaptivity to Hp-adaptivity

Test

$hp$-adaptivity

2022-09-20T05:40:15Z

Mitjajancic: Mitjajancic moved page $hp$-adaptivity to Hp-adaptivity

#REDIRECT [[Hp-adaptivity]]

HP-adaptivity

2022-09-20T05:39:55Z

Mitjajancic: Mitjajancic moved page Hp adaptivity to $hp$-adaptivity

Test

Hp adaptivity

2022-09-20T05:39:55Z

Mitjajancic: Mitjajancic moved page Hp adaptivity to $hp$-adaptivity

#REDIRECT [[$hp$-adaptivity]]

HP-adaptivity

2022-09-20T05:39:39Z

Mitjajancic: Created page with "Test"

Test

Plasticity

2021-10-08T10:23:13Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of plastic deformation by employing the von Mises plasticity model with isotropic hardening. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>. For a more detailed theory description, we strongly recommend reading the book.

= Introduction =
In physics and materials science, ''plasticity'', also known as ''plastic deformation'', is observed when the solid material undergoes permanent deformation in response to applied forces. There are different small-strain elastoplastic constitutive models of plasticity known. The most popular theories are the von Mises, Tresca, Mohr-Coulomb and Drucker-Prager models. The main difference between them is in the computation of yield condition - threshold value after which the deformation is irreversible.

This demonstration uses ''the von Mises model, with nonlinear isotropic hardening'', due to the simplicity of its computational implementation and appropriateness for real-world applications.

= Theory of plasticity =
In continuum mechanics, stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other. When additional external force $\bf f$ is applied to the body, the dynamics can be described by a second Newton's law
\begin{equation}
\ \nabla \cdot \boldsymbol{\sigma} + \bf{f} = 0,
\end{equation}
where $\boldsymbol \sigma $ is the stress tensor. Now contrary to the elastic response, where the relationship between the deformation $\bf u$ and stress $\boldsymbol \sigma$ is linear, the plastic response is a typical example of a non-linear relationship, because the stress tensor not only depends depends on the current deformation $\bf u$ but also on its history states.
Generally speaking, for such highly non-linear problems, iterative numerical methods are used to solve the problem. Therefore, the residuum vector $\bf r$ is introduced and above equation is rewritten
\begin{equation}
\ \bf r (\bf u_{n + 1}) = \nabla \boldsymbol \sigma (\bf u_{n + 1}) - \bf f = 0.
\end{equation}

TO BE FINISHED.
As previously stated, all of the theory is thoroughly explained in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto, mainly chapters 4-7.
=References=
<references/>

Plasticity

2021-10-08T10:22:33Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of plastic deformation by employing the von Mises plasticity model with isotropic hardening. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>. For a more detailed theory description, we strongly recommend reading the book.

= Introduction =
In physics and materials science, ''plasticity'', also known as ''plastic deformation'', is observed when the solid material undergoes permanent deformation in response to applied forces. There are different small-strain elastoplastic constitutive models of plasticity known. The most popular theories are the von Mises, Tresca, Mohr-Coulomb and Drucker-Prager models. The main difference between them is in the computation of yield condition - threshold value after which the deformation is irreversible.

This demonstration uses ''the von Mises model, with nonlinear isotropic hardening'', due to the simplicity of its computational implementation and appropriateness for real-world applications.

= Theory of plasticity =
In continuum mechanics, stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other. When additional external force $\bf f$ is applied to the body, the dynamics can be described by a second Newton's law
\begin{equation}
\ \nabla \cdot \boldsymbol{\sigma} + \bf{f} = 0,
\end{equation}
where $\boldsymbol \sigma $ is the stress tensor. Now contrary to the elastic response, where the relationship between the deformation $\bf u$ and stress $\boldsymbol \sigma$ is linear, the plastic response is a typical example of a non-linear relationship, because the stress tensor not only depends depends on the current deformation $\bf u$ but also on its history states.
Generally speaking, for such highly non-linear problems, iterative numerical methods are used to solve the problem. Therefore, the residuum vector $\bf r$ is introduced and above equation is rewritten
\begin{equation}
\ \bf r (\bf u_{n + 1}) = \nabla \boldsymbol \sigma (\bf u_{n + 1}) - \bf f = 0.
\end{equation}

=References=
<references/>

TO BE FINISHED.
As previously stated, all of the theory is thoroughly explained in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto, mainly chapters 4-7.

Plasticity

2021-10-08T10:20:56Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of plastic deformation by employing the von Mises plasticity model with isotropic hardening. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>. For a more detailed theory description, we strongly recommend reading the book.
<div class="mw-collapsible-content">
= Introduction =
In physics and materials science, ''plasticity'', also known as ''plastic deformation'', is observed when the solid material undergoes permanent deformation in response to applied forces. There are different small-strain elastoplastic constitutive models of plasticity known. The most popular theories are the von Mises, Tresca, Mohr-Coulomb and Drucker-Prager models. The main difference between them is in the computation of yield condition - threshold value after which the deformation is irreversible.

This demonstration uses ''the von Mises model, with nonlinear isotropic hardening'', due to the simplicity of its computational implementation and appropriateness for real-world applications.

= Theory of plasticity =
In continuum mechanics, stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other. When additional external force $\bf f$ is applied to the body, the dynamics can be described by a second Newton's law
\begin{equation}
\ \nabla \cdot \boldsymbol{\sigma} + \bf{f} = 0,
\end{equation}
where $\boldsymbol \sigma $ is the stress tensor. Now contrary to the elastic response, where the relationship between the deformation $\bf u$ and stress $\boldsymbol \sigma$ is linear, the plastic response is a typical example of a non-linear relationship, because the stress tensor not only depends depends on the current deformation $\bf u$ but also on its history states.
Generally speaking, for such highly non-linear problems, iterative numerical methods are used to solve the problem. Therefore, the residuum vector $\bf r$ is introduced and above equation is rewritten
\begin{equation}
\ \bf r (\bf u_{n + 1}) = \nabla \boldsymbol \sigma (\bf u_{n + 1}) - \bf f = 0.
\end{equation}

TO BE FINISHED.
As previously stated, all of the theory is thoroughly explained in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto, mainly chapters 4-7.

Plasticity

2021-10-08T10:20:07Z

Mitjajancic:

Plasticity

2021-10-08T10:01:07Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of '''plastic deformation by employing the von Mises plasticity model with isotropic hardening'''. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>. For a more detailed theory description, we strongly recommend reading the book.

= Introduction =
In physics and materials science, ''plasticity'', also known as plastic deformation, is observed when the solid material undergoes permanent deformation in response to applied forces. There are different small-strain elastoplastic constitutive models of plasticity known. The most popular theories are the von Mises, Tresca, Mohr-Coulomb and Drucker-Prager models. The main difference between them is in the computation of yield condition - threshold value after which the deformation transits from elastic to plastic range.

This demonstration uses ''the von Mises model, with nonlinear isotropic hardening'', due to its computational implementation simplicity and appropriateness for real-world examples.

= Theory of plasticity =
In continuum mechanics, stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other. When additional external force is applied to the body, the dynamics are described by a second Newton's law
\begin{equation}
\ \nabla \cdot \boldsymbol{\sigma} + \bf{f} = 0
\end{equation}

As previously stated, all of the theory is thoroughly explained in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>, mainly chapters 4-7.

Plasticity

2021-10-08T10:00:19Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of '''plastic deformation by employing the von Mises plasticity model with isotropic hardening'''. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>. For a more detailed theory description, we strongly recommend reading the book.

= Introduction =
In physics and materials science, ''plasticity'', also known as plastic deformation, is observed when the solid material undergoes permanent deformation in response to applied forces. There are different small-strain elastoplastic constitutive models of plasticity known. The most popular theories are the von Mises, Tresca, Mohr-Coulomb and Drucker-Prager models. The main difference between them is in the computation of yield condition - threshold value after which the deformation transits from elastic to plastic range.

This demonstration uses ''the von Mises model, with nonlinear isotropic hardening'', due to its computational implementation simplicity and appropriateness for real-world examples.

= Theory of plasticity =
In continuum mechanics, stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other. When additional external force is applied to the body, the dynamics are described by a second Newton's law
\begin{equation}
\ \nabla \cdot \bf{\sigma} + \bf{f} = 0
\end{equation}

As previously stated, all of the theory is thoroughly explained in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>, mainly chapters 4-7.

Plasticity

2021-10-08T10:00:09Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of '''plastic deformation by employing the von Mises plasticity model with isotropic hardening'''. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>. For a more detailed theory description, we strongly recommend reading the book.

= Introduction =
In physics and materials science, ''plasticity'', also known as plastic deformation, is observed when the solid material undergoes permanent deformation in response to applied forces. There are different small-strain elastoplastic constitutive models of plasticity known. The most popular theories are the von Mises, Tresca, Mohr-Coulomb and Drucker-Prager models. The main difference between them is in the computation of yield condition - threshold value after which the deformation transits from elastic to plastic range.

This demonstration uses ''the von Mises model, with nonlinear isotropic hardening'', due to its computational implementation simplicity and appropriateness for real-world examples.

= Theory of plasticity =
In continuum mechanics, stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other. When additional external force is applied to the body, the dynamics are described by a second Newton's law
\begin{equation}
\ \nabla \cdot \bf{sigma} + \bf{f} = 0
\end{equation}

As previously stated, all of the theory is thoroughly explained in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>, mainly chapters 4-7.

Plasticity

2021-10-08T09:28:42Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of '''plastic deformation with von Mises plasticity model and isotropic hardening'''. The majority of the theory is explained in detail in '''Computational Methods for Plasticity - Theory and Applications''' by EA de Souza Neto <ref> de Souza Neto, E.A., Peric, D. and Owen, D.R., 2011. Computational methods for plasticity: theory and applications. John Wiley & Sons</ref>.

Plasticity

2021-10-08T09:21:57Z

Mitjajancic:

Do you want to go back to [[Solid Mechanics]]?

On this page we conduct numerical studies of '''bending of a cantilever loaded at the end''', a common numerical benchmark in elastostatics. <ref> Augarde, Charles E. and Deeks, Andrew J.. "The use of Timoshenko's exact solution for a cantilever beam in adaptive analysis" , ''Finite Elements in Analysis and Design''. (2008), doi: [http://dx.doi.org/10.1016/j.finel.2008.01.010 10.1016/j.finel.2008.01.010] </ref>

Plasticity

2021-10-08T09:21:02Z

Mitjajancic: Created page with " On this page we showcase some basic examples from linear elasticity. To read more about the governing equations, refer to the Solid Mechanics page and examples therein, w..."

On this page we showcase some basic examples from linear elasticity. To read more about the governing equations, refer to the [[Solid Mechanics]] page
and examples therein, which are considered in more detail.

Medusa

2021-10-08T09:20:11Z

Mitjajancic: /* Discussions / Applications */

'''Welcome to the Medusa wiki. To visit the main website, go to [http://e6.ijs.si/medusa/ http://e6.ijs.si/medusa/].'''

In [http://e6.ijs.si/ParallelAndDistributedSystems/ Parallel and Distributed Systems Laboratory] we are working on a C++ library that is first and foremost focused on tools for solving Partial Differential Equations by meshless methods. The basic idea is to create generic codes for tools that are needed for solving not only PDEs but many other problems, e.g. Moving Least Squares approximation, $k$-d tree, domain generation engines, etc.
We call this open source meshless project [http://e6.ijs.si/medusa/ Medusa: Coordinate Free Meshless Method implementation (MM)].

Technical details about code and examples can be found on our [http://e6.ijs.si/medusa/docs/ documentation page] and [https://gitlab.com/e62Lab/medusa Gitlab repository]. [[File:C.png|100px||link=https://gitlab.com/e62Lab/medusa|alt=Alt text|code]] [[File:doxygen.png|100px|link=http://e6.ijs.si/medusa/docs/|alt=Alt text|Documentation page]]

This wiki site is meant for more relaxed discussions about general principles, possible and already implemented applications, preliminary analyses, etc.
Note, that there are many grammatical mistakes, typos, stupid sentences, etc. This wiki is meant for quick information exchange and therefore we do not invest a lot of energy into styling :).

== Documentation ==
* [https://gitlab.com/e62Lab/medusa Code on Gitlab]
* [[How to build | Installation and building]]
* [[Including this library in your project | Including this library in your project]]
* [[Testing | Running tests]]
* [http://e6.ijs.si/medusa/docs/ Technical documentation]
* [[Coding style | Coding style]]
* [[Wiki editing guide | Wiki editing and backup guide]]

== Examples ==
In this section we present exact examples. Each of the below solutions can be found also in in the repository under examples. More explanation about the physical background and solution procedure can be found in following sections.
* [[Philosophy of examples and how to run them]]
* [[Poisson's equation]]
* [[Heat equation]]
* [[Linear elasticity]]
* [[Complex-valued problems]]
* [[Coupled domains]]
* [[Parametric domains | Parametric domains – Curved surface with variable density]]
* [[NURBS domains | Domains modeled with non-uniform rational basis spline's (NURBS)]]
* [[Realistic 3D models]]
* [[Customization]]
* [[Ghost nodes]]
* [[Electromagnetic scattering]]
* [[Schrödinger equation]]
* [[Wave equation]]
* [[Cahn-Hilliard equation]]
* [[Non-Newtonian_fluid_example | Non-Newtonian fluid]]
* [[Meshless Lattice Boltzmann method]]

== Building blocks ==
Medusa is modular coordinate-free parallel implementation of a numerical framework designed, but not limited to, for solving PDEs. In this section we present main modules of the library that can be also used as a standalone tools.
* [[Positioning of computational nodes]]
* [[Relaxation of the nodal distribution]]
* [[Refinement of the nodal distribution]]
* [[k-d tree|''k''-d tree]] and other spatial search structures
* Solving [[Solving linear systems | linear systems]], [[Solving overdetermined systems | overdetermined]] and [[Solving underdetermined systems | underdetermined]]
* [[Weighted Least Squares (WLS)]]
* [[Computation of shape functions]]
* [[Meshless Local Strong Form Method (MLSM)]]
* [[Radial basis function-generated finite differences (RBF-FD)]]
* [[Ghost nodes (theory)]]
* [[Integrators for time stepping]]
* [[RBF Interpolation]]

== Discussions / Applications ==
This section is meant for general discussion about the physical background of the examples, the solution procedures, various applications, etc. Note, that code snippets presented in discussion might not reflect the actual state of Medusa.
* [[Convection Diffusion equation | Convection Diffusion equation]]
* [[Adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
** [[Plasticity]]
* [[Fluid Mechanics]]
** [[Lid driven cavity]]
** [[de Vahl Davis natural convection test]]
** [[Natural convection in 3D irregular domain]]
** [[Natural convection from heated cylinder]]
** [[Natural convection between concentric cylinders]]
** [[Non-Newtonian fluid]]
* [[Computational electromagnetics]]
** [[Triple dielectric step in 1D]]
** [[Scattering from an infinite cylinder]]
** [[Point source near an anisotropic lens]]
* Other applications
** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]
** [[Heart rate variability detection]]
** [[Bioheat equation]]

== Performance analyses ==
* [[Execution on Intel® Xeon Phi™ co-processor]]
* [[1D MLSM and FDM comparison]]
* [[:File:tech_report.pdf|Execution overheads due to clumsy types::technical report]] [[File:pdf-file.gif]]
* [[Solving sparse systems]]
* [[Eigen paralelization]]

== Last changes ==
<news unique=1 limit = 5>
*{{{timeanddate}}} :: {{{title}}}

</news>

== Miscellaneous ==
* FAQ - [[Frequently asked questions]].
* [[List of wiki contributors]]
* List of library contributors: [http://e6.ijs.si/medusa/about#about-contributors See the official website]

== References ==

For pre-prints check http://e6.ijs.si/ParallelAndDistributedSystems/publications/

* Slak J., Kosec G., On Generation of Node Distributions for Meshless PDE Discretizations, SIAM Journal on Scientific Computing, 41(5), A3202
* Slak J., Kosec G. Adaptive radial basis function-generated finite differences method for contact problems. International journal for numerical methods in engineering, ISSN 0029-5981 [http://www-e6.ijs.si/ParallelAndDistributedSystems/pdf/32230439.pdf manuscript]
* Slak J., Kosec G.; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain. Engineering analysis with boundary elements. 2018;11 ; [http://comms.ijs.si/~gkosec/data/papers/31107623.pdf manuscript]
* Depolli, M., Kosec, G., Assessment of differential evolution for multi-objective optimization in a natural convection problem solved by a local meshless method. Engineering optimization, 2017, vol. 49, no. 4, pp. 675-692 ;[http://comms.ijs.si/~gkosec/data/papers/29639719.pdf manuscript]
* Kosec G., A local numerical solution of a fluid-flow problem on an irregular domain. Advances in engineering software. 2016;7 ; [29512743] ; [http://comms.ijs.si/~gkosec/data/papers/29512743.pdf manuscript]
* Kosec G., Trobec R., Simulation of semiconductor devices with a local numerical approach. Engineering analysis with boundary elements. 2015;69-75; [27912487] ; [http://comms.ijs.si/~gkosec/data/papers/27912487.pdf manuscript]
* Kosec G., Šarler B., Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method. Engineering analysis with boundary elements. 2014;36-44; [http://comms.ijs.si/~gkosec/data/papers/3218939.pdf manuscript]
* Kosec G., Depolli M., Rashkovska A., Trobec R., Super linear speedup in a local parallel meshless solution of thermo-fluid problem. Computers & Structures. 2014;133:30-38; [http://comms.ijs.si/~gkosec/data/papers/27339815.pdf manuscript]
* Kosec G., Zinterhof P., Local strong form meshless method on multiple Graphics Processing Units. Computer modeling in engineering & sciences. 2013;91:377-396; [http://comms.ijs.si/~gkosec/data/papers/26785063.pdf manuscript]
* Kosec G., Šarler B., H-adaptive local radial basis function collocation meshless method. Computers, materials & continua. 2011;26:227-253; [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerBurgers.pdf manuscript]
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321; [http://comms.ijs.si/~gkosec/data/papers/EABE2499.pdf manuscript]
* Kosec G, Zaloznik M, Sarler B, Combeau H. A Meshless Approach Towards Solution of Macrosegregation Phenomena. CMC: Computers, Materials, & Continua. 2011;580:1-27 [http://comms.ijs.si/~gkosec/data/papers/KosecZaloznikSarlerCombeauSegregation.pdf manuscript]
* Kosec G, Sarler B. Solution of thermo-fluid problems by collocation with local pressure correction. International Journal of Numerical Methods for Heat & Fluid Flow. 2008;18:868-82 [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerNS2008.pdf manuscript]
* Trobec R., Kosec G., Parallel Scientific Computing, ISBN: 978-3-319-17072-5 (Print) 978-3-319-17073-2.
* Slak, J., Kosec, G.. Detection of heart rate variability from a wearable differential ECG device., MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938, pp 450-455.
* Kolman, M., Kosec, G. Correlation between attenuation of 20 GHz satellite communication link and liquid water content in the atmosphere. MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938. pp. 308-313.

==Related pages==
* http://e6.ijs.si/ParallelAndDistributedSystems/products/medusa/

Coding style

2020-11-30T13:06:11Z

Mitjajancic: /* Using your IDE auto formatter */

This is a brief description of our coding style, roughtly following the [https://google.github.io/styleguide/cppguide.html Google C++ Style Guide].

You can also download the CLion settings: [[File:Settings.jar|thumb]], which you can import in <code>File/Import settings</code>.

You can also download only our code style settings to import in CLion: [[:File:mm.xml]]. It can be imported under <code>Settings/Editor/Code Style/Scheme/Import scheme</code>.

=General=
We use a 100 characters line width limit in C++ code.

=Indentation=

Indent using spaces. Indentation width is 4 spaces.
<syntaxhighlight lang="c++">
for (int i = 0; i < 10; i++) {
cout <

=Naming convention=

Constants - <syntaxhighlight lang="c++" inline> UPPER_CASE </syntaxhighlight>

classes - <syntaxhighlight lang="c++" inline>PascalCase</syntaxhighlight>

methods - <syntaxhighlight lang="c++" inline>camelCase</syntaxhighlight>

variables and stand-alone functions - <syntaxhighlight lang="c++" inline> underscore_separated </syntaxhighlight>

typedefs - lowercase underscore separated, usually one word with trailing <syntaxhighlight lang="c++" inline> _t </syntaxhighlight>, eg. <syntaxhighlight lang="c++" inline> vec_t </syntaxhighlight>.

namespaces - one lowercase word, maybe shortened, eg. <syntaxhighlight lang="c++" inline> op </syntaxhighlight>. For internal implementation details use <syntaxhighlight lang="c++" inline> op_internal </syntaxhighlight>. Namespace closing brace should contain comment `// namespace your_name`. There is no indentation within namespaces.

All standard abbreviations like [[Weighted Least Squares (WLS)]] or [https://en.wikipedia.org/wiki/Finite_element_method Finite Elements Method (FEM)] -
<syntaxhighlight lang="c++" inline>UPPER_CASE</syntaxhighlight>

Floating point and integer literals use small suffixes, eg. <syntaxhighlight lang="c++" inline> 0.0f, -1e8l, 45ull</syntaxhighlight>
.

<syntaxhighlight lang="c++">
#define MAX_VALUE 500
class MyClass {
public:
MyClass(first_var, second_var) {
...
}
int getSize();
}
</syntaxhighlight>

=Comments=

Comments are good. Use them to explain your code. Comments should have a space between the last
slash and the start of text. Inline comments should have at least two spaces between end of code
and start of the comment. Comments before functions are used for documentation. Always use full
sentences and end them with a punctuation mark.

<syntaxhighlight lang="c++">
/**
* This function changes the world.
* @param skynet If `true`, the world ends, otherwise nothing happens.
*/
double change_the_world(bool skynet) {
if (skynet) {
return 0.0; // Brace for the end of the world.
}
...
}
</syntaxhighlight>

Use doxygen comments to generate [http://e6.ijs.si/medusa/docs/html/ documentation].

=Headers =

All headers must contain a header guard of form <syntaxhighlight lang="c++" inline>PATH_TO_FILENAME_HPP_</syntaxhighlight> as enforced by the linter.

Includes in header guards are separated in two groups, with intra-project includes on top and other includes on the bottom.
The groups are separated by a blank line and includes are kept sorted within a group.

=Misc=

Avoid trailing whitespace. Curly opening brackets <syntaxhighlight lang="c++" inline> { </syntaxhighlight> should be inline with <syntaxhighlight lang="c++" inline> for </syntaxhighlight> loops, function
definitions, class names, separated with a space. Outermost binary operators should have spaces
around them.

<syntaxhighlight lang="c++">
int sumMe(int var) { // yes
if (var == 1)
{ // no
return 1;
}
return 0;
}
</syntaxhighlight>
For null pointer we use <syntaxhighlight lang="c++" inline> nullptr </syntaxhighlight> instead of <syntaxhighlight lang="c++" inline> NULL </syntaxhighlight> macro.

=Using your IDE auto formatter=
Most available IDEs support source code formatting using [https://clang.llvm.org/docs/ClangFormat.html clang-format]. We recommend using coding style based on Google, with minor modifications concerning indentation and pointer location.

<syntaxhighlight>
# Use the Google style in this project.
BasedOnStyle: Google

# Custom adjustments.
AccessModifierOffset: -2
IndentWidth: 4
ColumnLimit: 100

# Some folks prefer to write "int& foo" while others prefer "int &foo". The
# Google Style Guide only asks for consistency within a project, we chose
# "int& foo" for this project:
DerivePointerAlignment: false
PointerAlignment: Left

# Do not sort includes.
SortIncludes: false
</syntaxhighlight>

More information about auto formatting in:
* [https://code.visualstudio.com/docs/cpp/cpp-ide Visual Studio Code],
* [https://www.jetbrains.com/help/clion/clangformat-as-alternative-formatter.html Jetbrains Clion]

Coding style

2020-11-30T13:00:08Z

Mitjajancic: Updated .clang-format.

This is a brief description of our coding style, roughtly following the [https://google.github.io/styleguide/cppguide.html Google C++ Style Guide].

You can also download the CLion settings: [[File:Settings.jar|thumb]], which you can import in <code>File/Import settings</code>.

You can also download only our code style settings to import in CLion: [[:File:mm.xml]]. It can be imported under <code>Settings/Editor/Code Style/Scheme/Import scheme</code>.

=General=
We use a 100 characters line width limit in C++ code.

=Indentation=

Indent using spaces. Indentation width is 4 spaces.
<syntaxhighlight lang="c++">
for (int i = 0; i < 10; i++) {
cout <

=Naming convention=

Constants - <syntaxhighlight lang="c++" inline> UPPER_CASE </syntaxhighlight>

classes - <syntaxhighlight lang="c++" inline>PascalCase</syntaxhighlight>

methods - <syntaxhighlight lang="c++" inline>camelCase</syntaxhighlight>

variables and stand-alone functions - <syntaxhighlight lang="c++" inline> underscore_separated </syntaxhighlight>

typedefs - lowercase underscore separated, usually one word with trailing <syntaxhighlight lang="c++" inline> _t </syntaxhighlight>, eg. <syntaxhighlight lang="c++" inline> vec_t </syntaxhighlight>.

namespaces - one lowercase word, maybe shortened, eg. <syntaxhighlight lang="c++" inline> op </syntaxhighlight>. For internal implementation details use <syntaxhighlight lang="c++" inline> op_internal </syntaxhighlight>. Namespace closing brace should contain comment `// namespace your_name`. There is no indentation within namespaces.

All standard abbreviations like [[Weighted Least Squares (WLS)]] or [https://en.wikipedia.org/wiki/Finite_element_method Finite Elements Method (FEM)] -
<syntaxhighlight lang="c++" inline>UPPER_CASE</syntaxhighlight>

Floating point and integer literals use small suffixes, eg. <syntaxhighlight lang="c++" inline> 0.0f, -1e8l, 45ull</syntaxhighlight>
.

<syntaxhighlight lang="c++">
#define MAX_VALUE 500
class MyClass {
public:
MyClass(first_var, second_var) {
...
}
int getSize();
}
</syntaxhighlight>

=Comments=

Comments are good. Use them to explain your code. Comments should have a space between the last
slash and the start of text. Inline comments should have at least two spaces between end of code
and start of the comment. Comments before functions are used for documentation. Always use full
sentences and end them with a punctuation mark.

<syntaxhighlight lang="c++">
/**
* This function changes the world.
* @param skynet If `true`, the world ends, otherwise nothing happens.
*/
double change_the_world(bool skynet) {
if (skynet) {
return 0.0; // Brace for the end of the world.
}
...
}
</syntaxhighlight>

Use doxygen comments to generate [http://e6.ijs.si/medusa/docs/html/ documentation].

=Headers =

All headers must contain a header guard of form <syntaxhighlight lang="c++" inline>PATH_TO_FILENAME_HPP_</syntaxhighlight> as enforced by the linter.

Includes in header guards are separated in two groups, with intra-project includes on top and other includes on the bottom.
The groups are separated by a blank line and includes are kept sorted within a group.

=Misc=

Avoid trailing whitespace. Curly opening brackets <syntaxhighlight lang="c++" inline> { </syntaxhighlight> should be inline with <syntaxhighlight lang="c++" inline> for </syntaxhighlight> loops, function
definitions, class names, separated with a space. Outermost binary operators should have spaces
around them.

<syntaxhighlight lang="c++">
int sumMe(int var) { // yes
if (var == 1)
{ // no
return 1;
}
return 0;
}
</syntaxhighlight>
For null pointer we use <syntaxhighlight lang="c++" inline> nullptr </syntaxhighlight> instead of <syntaxhighlight lang="c++" inline> NULL </syntaxhighlight> macro.

=Using your IDE auto formatter=
Most available IDEs support source code formatting using [https://clang.llvm.org/docs/ClangFormat.html clang-format]. We recommend using coding style based on Google, with minor modifications concerning indentation and pointer location.

<syntaxhighlight>
# Use the Google style in this project.
BasedOnStyle: Google

# Custom adjustments.
AccessModifierOffset: -2
IndentWidth: 4
ColumnLimit: 100

# Some folks prefer to write "int& foo" while others prefer "int &foo". The
# Google Style Guide only asks for consistency within a project, we chose
# "int& foo" for this project:
DerivePointerAlignment: false
PointerAlignment: Left

# Do not sort includes.
SortIncludes:false
</syntaxhighlight>

More information about auto formatting in:
* [https://code.visualstudio.com/docs/cpp/cpp-ide Visual Studio Code],
* [https://www.jetbrains.com/help/clion/clangformat-as-alternative-formatter.html Jetbrains Clion]