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

Wave equation

2019-05-08T11:14:28Z

JMocnik: Terminology corrected.

==2D wave equation with Dirichlet boundary conditions ==
We are solving a problem analogous to a clamped circular membrane being forced to oscillate by fixing its centre to a rigid oscillating rod. The domain of calculation is for this reason an annulus bounded by the point of contact with the circular source in the centre and the clamped edge of the circular membrane.

[[File:wave_2d_domain_geo.jpg|300px]]

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u(r,t=0)&=0 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in accordance with the desired density function. Once the domain is constructed and discretized, the support (neighborhood) of $n$ nodes is found for each node.

<syntaxhighlight lang="c++" line>

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

The density function "fill_density" takes as argument the coordinates inside the domain and returns the desired distance between neighboring nodes for that location.
Now the domain preparation is complete, as domain is constructed and populated with nodes with support relations established. We continue by constructing the approximation engine.
<syntaxhighlight lang="c++" line>
WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> approx(m - 1 , sigma);
</syntaxhighlight>

As we are solving a time propagation problem implicitly, we will have to solve a system of equations for every time step. This is equivalent to solving a matrix equation. For this reason a space matrix is constructed. Next the wave equation is turned into code. Medusa encodes this information inside the just constructed matrix $M$. $E1$ and $E0$ represent the known solutions on previous two time steps.

<syntaxhighlight lang="c++" line>
SparseMatrix<double> M(domain_size, domain_size);
// Set equation on interior
for (int i : interior) {
-(v*v*dt*dt) * op.lap(i) + op.value(i)= 2*E1(i)-E0(i);
}
</syntaxhighlight>
Same is done also for the boundary nodes. Since dirichlet boundary conditions are employed the values of the Right Hand Side ($rhs$) of matrix equation are just set to prescribed values ($0$ on the outer boundary and $f(t)$ on inner boundary).
<syntaxhighlight lang="c++" line>
// Set boundary conditions
for (int i : boundary) {
if (discretization.types()[i] == CENTRE) {
op.value(i) = A * std::sin(omega*dt*2); // step
} else {
op.value(i) = 0.0;
}
}
</syntaxhighlight>
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BiCGStab algorithm. $E0$ and $E1$ are updated and the source is set to new value.

<syntaxhighlight lang="c++" line>
// Solve matrix system
VectorXd E2 = solver.solve(rhs);

// Update previous states
E0 = E1;
E1 = E2;

// Update rhs on interior
for (int i : interior) {
rhs(i) = 2*E1(i)-E0(i);
}

// Update rhs on boundary
for (int i : (discretization.types() == CENTRE)) {
rhs(i) = A * std::sin(omega*dt*tt); // step
}
</syntaxhighlight>
The animation of solution is presented below. The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/tree/dev/examples/wave_equation wave_equation_2D.cpp]. This example is covered in more detail in linked document [[:File:Meshless solution to wave equation.pdf]].

[[File:EM_wave_2d_dt_e-6_dens40_omega128_v8 2X speed.mp4|400px]]

Wave equation

2019-05-08T11:06:56Z

JMocnik: Example desctripction completed

==2D wave equation with Dirichlet boundary conditions ==
We are solving a problem analogous to a clamped circular membrane being forced to oscillate by fixing its centre to a rigid oscillating rod. The domain of calculation is for this reason an annulus bounded by the point of contact with the circular source in the centre and the clamped edge of the circular membrane.

[[File:wave_2d_domain_geo.jpg|300px]]

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u(r,t=0)&=0 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in accordance with the desired density function. Once the domain is constructed and discretized, the support (neighborhood) of $n$ nodes is found for each node.

<syntaxhighlight lang="c++" line>

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

The density function "fill_density" takes as argument the coordinates inside the domain and returns the desired distance between neighboring nodes for that location.
Now the domain preparation is complete, as domain is constructed and populated with nodes with support relations established. We continue by constructing the approximation engine.
<syntaxhighlight lang="c++" line>
WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> approx(m - 1 , sigma);
</syntaxhighlight>

As we are solving a time propagation problem implicitly, we will have to solve a system of equations for every time step. This is equivalent to solving a matrix equation. For this reason a space matrix is constructed. Next the wave equation is turned into code. Medusa encodes this information inside the just constructed matrix $M$. $E1$ and $E0$ represent the known solutions on previous two time steps.

<syntaxhighlight lang="c++" line>
SparseMatrix<double> M(domain_size, domain_size);
// Set equation on interior
for (int i : interior) {
-(v*v*dt*dt) * op.lap(i) + op.value(i)= 2*E1(i)-E0(i);
}
</syntaxhighlight>
Same is done also for the boundary nodes. Since dirichlet boundary conditions are employed the values of the Right Hand Side ($rhs$) of matrix equation are just set to prescribed values ($0$ on the exterior and $f(t)$ on interior).
<syntaxhighlight lang="c++" line>
// Set boundary conditions
for (int i : boundary) {
if (discretization.types()[i] == CENTRE) {
op.value(i) = A * std::sin(omega*dt*2); // step
} else {
op.value(i) = 0.0;
}
}
</syntaxhighlight>
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BiCGStab algorithm. $E0$ and $E1$ are updated and the source is set to new value.

<syntaxhighlight lang="c++" line>
// Solve matrix system
VectorXd E2 = solver.solve(rhs);

// Update previous states
E0 = E1;
E1 = E2;

// Update rhs on interior
for (int i : interior) {
rhs(i) = 2*E1(i)-E0(i);
}

// Update rhs on boundary
for (int i : (discretization.types() == CENTRE)) {
rhs(i) = A * std::sin(omega*dt*tt); // step
}
</syntaxhighlight>
The animation of solution is presented below. The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/tree/dev/examples/wave_equation wave_equation_2D.cpp]. This example is covered in more detail in linked document [[:File:Meshless solution to wave equation.pdf]].

[[File:EM_wave_2d_dt_e-6_dens40_omega128_v8 2X speed.mp4|400px]]

Wave equation

2019-05-08T09:17:34Z

JMocnik:

==2D wave equation with Dirichlet boundary conditions ==
We are solving a problem analogous to a clamped circular membrane being forced to oscillate by fixing its centre to a rigid oscillating rod. The domain of calculation is for this reason an annulus bounded by the point of contact with the circular source in the centre and the clamped edge of the circular membrane.

[[File:wave_2d_domain_geo.jpg|300px]]

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u(r,t=0)&=0 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in accordance with the desired density function. Once the domain is constructed and discretized, the support (neighborhood) of $n$ nodes is found for each node.

<syntaxhighlight lang="c++" line>

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

The density function "fill_density" takes as argument the coordinates inside the domain and returns the desired distance between neighboring nodes for that location.
Now the domain preparation is complete, as domain is constructed and populated with nodes with support relations established. We continue by constructing the approximation engine.
<syntaxhighlight lang="c++" line>
WLS<Monomials<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> approx(m - 1 , sigma);
</syntaxhighlight>

As we are solving a time propagation problem implicitly, we will have to solve a system of equations for every time step. This is equivalent to solving a matrix equation. For this reason a space matrix is constructed. Next the wave equation is turned into code. Medusa encodes this information inside the just constructed matrix $M$. $E1$ and $E0$ represent the known solutions on previous two time steps.

<syntaxhighlight lang="c++" line>
SparseMatrix<double> M(domain_size, domain_size);
// Set equation on interior
for (int i : interior) {
M.coeffRef(i, i) = 1;
-(v*v*dt*dt) * op.lap(i) = 2*E1(i)-E0(i); // laplace in interior
}
</syntaxhighlight>
Same is done also for the boundary nodes. Since dirichlet boundary conditions are employed the values of the Right Hand Side ($rhs$) of matrix equation are just set to prescribed values ($0$ on the exterior and $f(t)$ on interior).
<syntaxhighlight lang="c++" line>
// Set boundary conditions
for (int i : boundary) {
M.coeffRef(i, i) = 1; // fixed
if (discretization.types()[i] == CENTRE) {
rhs(i) = A * std::sin(omega*dt*2); // step
} else {
rhs(i) = 0.0;
}
}
</syntaxhighlight>
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BiCGStab algorithm. $E0$ and $E1$ are updated and the source is set to new value.

<syntaxhighlight lang="c++" line>
// Solve matrix system
VectorXd E2 = solver.solve(rhs);

// Update previous states
E0 = E1;
E1 = E2;

// Update rhs on interior
for (int j = 0; j < interior.size(); ++j) {
int i = interior[j];
rhs(i) = 2*E1(i)-E0(i);
}

// Update rhs on boundary
for (int i : (discretization.types() == CENTRE)) {
rhs(i) = A * std::sin(omega*dt*tt); // step
}
</syntaxhighlight>
The animation of solution is presented below. The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/tree/dev/examples/wave_equation wave_equation_2D.cpp]. This example is covered in more detail in linked document: [[:File:Meshless solution to wave equation.pdf]].

[[File:EM_wave_2d_dt_e-6_dens40_omega128_v8 2X speed.mp4|400px]].

.

Wave equation

2019-05-08T07:13:01Z

JMocnik:

Wave equation

2019-05-07T17:10:03Z

JMocnik:

==2D wave equation with dirtchlet boundary conditions ==

Consider the time dependent solution to 2D wave equation on annulus shaped domain

<math>
\begin{align*}
\frac{ \partial^2 u}{\partial t^2} &= c^2 \nabla^2 u &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega_O,\\
u &= f(t) &&\text{on } \partial \Omega_I,
\end{align*}
</math>

where $\partial\Omega_I$ denotes the inner and $\partial\Omega_O$ the outer boundary of the domain. Through the boundary condition on the inner boundary the source is introduced to the problem as a function of time

<math>
\begin{align*}
f(t)= u_o \sin \omega_o t.
\end{align*}
</math>

First the domain is constructed by subtracting a smaller circle domain from a larger one. Boundaries of the domain are populated in the same step.

<syntaxhighlight lang="c++" line>
// // identifier to be added to nodes on the inner boundary
int CENTRE = -10;

// build circle domain
BallShape<Vec2d> domain({0, 0}, outer_radius);
auto discretization = domain.discretizeBoundaryWithStep(dx);

// build source domain
BallShape<Vec2d> empty({0, 0}, inner_radius);
auto discretization_empty = empty.discretizeBoundaryWithStep(dx, CENTRE);

// substract the source domain
discretization -= discretization_empty;
</syntaxhighlight>

Next the domain is populated with nodes in acordance with the desired density function.

<syntaxhighlight lang="c++" line>
// Lambda function for setting the density of nodes
auto fill_density = [=](const Vec2d& p) {
double r = p.norm();
double default_value = dx;
double dens = default_value;
double r1 = 15*inner_radius;
double r2 = 0.8*outer_radius;
if (r < r1) dens = linear(inner_radius, 0.8*default_value, r1, default_value, r );
if (r > r2) dens = linear(r2, default_value, outer_radius, 0.8* default_value, r);
return dens;
};

GeneralFill<Vec2d> fill;
fill.seed(fill_seed);
discretization.fill(fill, fill_density);

// find support
FindClosest find_support(n);
discretization.findSupport(find_support);
</syntaxhighlight>

Medusa

2018-08-29T18:31:46Z

JMocnik: Fixes tipo in Meshless

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 '''Medusa: Coordinate Free Meshless Method implementation (MM).

Technical details about code, examples, and can be found on our [http://e6.ijs.si/medusa/docs/ documentation page] and [https://gitlab.com/e62Lab/medusa the code]. [[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]]
* [[Testing | Running unit tests]]
* [http://e6.ijs.si/medusa/docs/ Technical documentation]
* [[Coding style | Coding style]]
* [[Including this library in your project | Including this library in your project]]
* [[Wiki editing guide | Wiki editing and backup guide]]

== Examples (work in progress) ==
* [[Philosophy of examples and how to run them]]
* [[Poisson's equation]]
* [[Linear elasticity]]
* [[Fluid mechanics]]
* [[Wave equation]]

== Building blocks ==
* [[Moving Least Squares (MLS)]]
* [[k-d tree|''k''-d tree]] TODO: UPDATE
* [[Positioning of computational nodes]]
* [[Meshless Local Strong Form Method (MLSM)]]
* [[Computation of shape functions]]
* [[Integrators for time stepping]]

== Applications / Examples ==
* Basic PDE solutions
** [[Heat Poisson's equation | Convection Diffusion equation]]
** [[Quantum Mechanics | Schrödinger equation]]
**[[Wave equation application]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[FWO case]]
* [[Fluid Mechanics]]
** [[Lid driven cavity]]
** [[de Vahl Davis natural convection test]]
** [[Natural convection from heated cylinder]]
** [[Natural convection between concentric cylinders]]
* Other applications
** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]
** [[Heart rate variability detection]]

== 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>

== MISC ==
* Also see FAQ -[[Frequently asked questions]].
* [[List of contributors]]

== Contributors ==
* TODO

== References ==
* 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://www-e6.ijs.si/ParallelAndDistributedSystems/#!NumericalMethods
* http://www-e6.ijs.si/ParallelAndDistributedSystems/#!utils
* http://www-e6.ijs.si/ParallelAndDistributedSystems/#!NUMA

File:EM wave 2d dt e-6 dens40 omega128 v8 2X speed.mp4

2018-04-20T09:13:30Z