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

Philosophy of examples and how to run them

2022-03-02T12:36:23Z

MihaR: Renamed the example file due to changes in the code.

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

The idea of our examples is to present the functionality of the library in an accessible and incremental manner.
The examples present different solutions to problems ranging in difficulty from the simple [[Poisson's equation]] to more specific examples, such as [[Fluid mechanics]].
Each examples consists of a name, e.g. <code>poisson_dirichlet_2D</code> and its corresponding
* source code, e.g. <code>poisson_dirichlet_2D.cpp</code>
* plot script, e.g. <code>poisson_dirichlet_2D.m</code>
Additional appropriately named input files might be present for more complicated examples and similarly named output files might be produced.
The examples and their plot scripts are intended to work out of the box, so the reader can simply run the examples on their machine.
The procedure for running examples is simple and the same for all examples. The accompanying plot scripts can be run with Matlab or octave and
they produce the figures the reader sees in the Examples category of this wiki.

Let us now illustrate the process of running an example <code>poisson_dirichlet_2D</code>.

First the example needs to be built (compiled). The building is done in the same <code>build/</code> directory
where medusa was build. To build and run the example, simply execute
<syntaxhighlight lang="bash">
make poisson_dirichlet_2D_run # builds and runs your example
</syntaxhighlight>
to compile and run your example. On first compilation, the whole library might compile as well, if not compiled previously, so it might take longer than usual.
Compiled binary used to run the example is placed next to the example source file. If only <code>make poisson_dirichlet_2D</code> (without _run) the
example is built, but not run.

To build all examples, you can just run
<syntaxhighlight lang="bash">
make examples # or "make examples_run" to also run them
</syntaxhighlight>

To run the built example again, go back to the initial folder where the source code for your example was, and run
<syntaxhighlight lang="bash">
./poisson_dirichlet_2D
</syntaxhighlight>
The output will be saved in the same folder as the binary and the source file.

After, just run the corresponding plot script <code>example_name.m</code> and you should see a figure, usually showing the solution.

'''Standalone compilation'''

The examples can then be compiled manually using e.g. <code>g++</code>, provided that
appropriate include and library paths are set and the <code>medusa_standalone</code> library is specified
for linking.

Example command to be run from <code>poisson_equation</code> directory:
<syntaxhighlight lang="bash">
g++ -o poisson_dirichlet_2D poisson_dirichlet_2D.cpp -I ../../include -Wall -O3 -L ../../bin -lmedusa_standalone
</syntaxhighlight>

The compiled example can be run simply as
<syntaxhighlight lang="bash">
./poisson_dirichlet_2D
</syntaxhighlight>
and plotted using the accompanying plot script.

Some more feature specific examples are found in the [http://e6.ijs.si/docs/ technical documentation] and in the unit tests.

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

How to build

2021-08-17T10:47:44Z

MihaR: Updated run script flags

=Installation=
To make this work from plain Ubuntu installation, run
<syntaxhighlight lang="bash">
sudo apt-get install git g++ python3 cmake libhdf5-serial-dev doxygen graphviz
git clone https://gitlab.com/e62Lab/medusa.git --branch master --single-branch
cd medusa
./run_tests.py -t
</syntaxhighlight>
which installs dependencies, clones the repository, goes into the root folder of
the repository and runs tests. This will build and run all tests. If this
works, you are ready to go! Otherwise install any missing packages and if it
still fails, raise an issue!

For instructions on how to use this library in you project, see
[[Including this library in your project]].

=Building=

List of dependencies:

* Build tools, like <code> cmake >= 2.8.12</code>, <code> g++ >= 4.8</code>, <code>make</code>, <code>python3</code>
* [https://www.hdfgroup.org/ HDF5 library] for IO
* <code>doxygen >= 1.8.8 </code> and Graphviz for generating the documentation

Out of source builds are preferred. Run
<syntaxhighlight lang="bash">
mkdir -p build
cd build
cmake ..
make
</syntaxhighlight>
Note that you only have to run <code>cmake</code> once, after that only <code>make</code> is sufficient.

Binaries are placed into <code>bin/</code> folder. Tests can be run all at once via <code>make medusa_run_tests</code>
or individually via e. g. <code>make operators_run_tests</code>.

== Linker errors ==

When trying out different classes, you might come along linker errors such as

<code>
Scanning dependencies of target cantilever_beam
[100%] Building CXX object examples/linear_elasticity/CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o
[100%] Linking CXX executable ../../../examples/linear_elasticity/cantilever_beam
/usr/bin/ld: CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o: in function `main':
cantilever_beam.cpp:(.text.startup+0x162): undefined reference to `void mm::FindBalancedSupport::operator()<mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> > >(mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&) const'
collect2: error: ld returned 1 exit status
</code>

This is expected and is the result of some optimizations of compilation time. The medusa library can actually be included in two ways: as
<code>#include <medusa/Medusa_fwd.hpp></code> or <code>#include <medusa/Medusa.hpp></code>. The first version contains the declarations of all symbols, but not all the definitions. Some of the more commonly used template instantiations are included, but by far not all. Using a template instantiation that is not precompiled will cause your program to compile fine, but will fail to link, due to the missing definitions. In this case you have a few options: include the <i>full</i> Medusa library (the header without the <code>_fwd</code>) and it should just work, but you will have to wait a bit longer for it to compile. Include only the missing header (in the case above <code>medusa/bits/domains/FindBalancedSupport.hpp</code>) and pay for whay you use. Or, add your instantiation among the already precompiled instantiations (located in <code>.cpp</code> files, such as e.g. [https://gitlab.com/e62Lab/medusa/-/blob/dev/src/domains/DomainDiscretization.cpp this one]).

== Building on macOS ==
This method was tested on macOS Mojave 10.14.2.

First install Xcode via App Store and then Xcode Command Line Tools with
<syntaxhighlight lang="bash">
xcode-select --install
</syntaxhighlight>

After that, install all dependencies from homebrew
<syntaxhighlight lang="bash">
brew install cmake hdf5 boost python doxygen graphviz
</syntaxhighlight>

Now you can clone and build the project using the following commands
<syntaxhighlight lang="bash">
git clone https://gitlab.com/e62Lab/medusa.git
cd medusa
mkdir build
cd build
cmake ..
cd ..
python3 run_tests.py -t
</syntaxhighlight>

This will also run all tests. If it works, you are ready to go! Otherwise install any missing packages and if it still fails, raise an issue!

==HDF5==
In order to use HDF5 IO you need the [https://www.hdfgroup.org/ HDF5 library].
You can install it easily using the command <code>sudo apt-get install libhdf5-dev</code>
or <code>sudo pacman -S hdf5</code>.

Ubuntu places (at least on older versions) hdf5 headers and libraries in a weird folder
<code>/usr/{lib, include}/x86_64-linux-gnu/hdf5/serial/</code>.

If you get an error like <code>HDF5 sample failed to compile. See errors above.</code> during <code>cmake</code> execution
then the sample hdf test file located in <code>test/test_hdf_compile.cpp</code> failed to compile. Perhaps it is good to make this file compile first,
before tackling the whole project.

If you get an error similar to <code>fatal error: hdf5.h: No such file or directory</code>,
then your compiler cannot see the HDF5 header files. Some distributions, notably (older) Ubuntu, place them into nonstandard folders
<code>/usr/include/hdf5/serial/</code> or <code>/usr/include/x86_64-linux-gnu/hdf5/serial/</code>.
Check these two folders or check your distributions hdf package for the locations of these files.
After you have determined the location, add that directory to the include directories,
using -I flag or in <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
include_directories(/usr/include/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
If you wish to fix this problem permanently, you can create soft links to the headers in your <code>/usr/include</code> directory,
by typing
<syntaxhighlight lang="bash">
sudo ln -s /usr/include/hdf5/serial/* /usr/include
</syntaxhighlight>
After this, there should be no compile time errors. If there are, please raise an issue.

If you get error similar to <code>-lhdf5 not found</code> and you have hdf5 installed,
you might have to link the libraries into a discoverable place, like <code>/usr/lib/</code>
or add the above directory to the linker path. Similarly to above, check the <code>/usr/lib/x86_64-linux-gnu/hdf5/serial/</code>
directory and look for file <code>libhdf5.a</code>. When found,
specify the location using -L flag or <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
link_directories(/usr/lib/x86_64-linux-gnu/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
or fix the problem permanently by soft-linking
<syntaxhighlight lang="bash">
sudo ln -s /usr/lib/x86_64-linux-gnu/hdf5/serial/* /usr/lib
</syntaxhighlight>

==Linear Algebra==
We use [http://eigen.tuxfamily.org/ Eigen] as our matrix library. See
[http://eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html here] for use
reference and documentation. For a quick transition from Matlab see
[http://eigen.tuxfamily.org/dox/AsciiQuickReference.txt here].

== Using Intel Math Kernel Library (MKL) ==
Install [https://software.intel.com/en-us/mkl Intel MKL] and take not of installation directories.

Proper include and link directories need to be set to use the Intel MKL.
<syntaxhighlight lang="cmake">
include_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/include) # change these to your installation path
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/lib/intel64)
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/lib/intel64)
</syntaxhighlight>

Eigen has great support for MKL. You can see more detailed instructions [https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html on their website].
To use MKL for math operations, define <code>EIGEN_USE_MKL_VML</code> when compiling. You must also link
the appropriate libraries and define <code>MKL_LP64</code> for use on a 64bit system.
With parallelism enabled, the configuration looks as follows

<syntaxhighlight lang="cmake">
target_compile_options(my_target PRIVATE "-fopenmp")
target_compile_definitions(my_target PUBLIC EIGEN_USE_MKL_VML MKL_LP64)
target_link_libraries(my_target medusa mkl_intel_lp64 mkl_intel_thread mkl_core pthread iomp5)
</syntaxhighlight>

If you have Intel Parallel Studio installed this also enables you to use the Pardiso paralle direct sparse solver through its [https://eigen.tuxfamily.org/dox/group__PardisoSupport__Module.html Eigen interface].

== Using Intel C/C++ Compiler ==

In order to use Intel's compiler when compiling Medusa, you have several standard optionas for <code>cmake</code>.
Make sure compilers and installed and in your <code>PATH</code> by running <code>which icc</code>, which
should return something like <code>/opt/intel/bin/icc</code>.

You can define the compiler when *first* calling cmake like so
<syntaxhighlight lang="bash">
cmake .. -DCMAKE_C_COMPILER=$(which icc) -DCMAKE_CXX_COMPILER=$(which icpc)
</syntaxhighlight>
If this is not your first call, remove the <code>build</code> directory and start anew.

You can also set the <code>CXX</code> and <code>CC</code> bash variables. Before calling
<code> cmake </code> for the first time you have to export the following:
<syntaxhighlight lang="bash">
export CXX="icpc"
export CC="icc"
</syntaxhighlight>

Fluid mechanics

2021-05-13T09:34:14Z

MihaR: Redirected page to Fluid Mechanics

#REDIRECT [[Fluid Mechanics]]

Medusa

2020-10-21T11:33:32Z

MihaR: /* Examples */ Added link to non-Newtonian example

'''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)]]
* [[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.
* Basic PDE solutions
** [[Convection Diffusion equation | Convection Diffusion equation]]
**[[Wave equation application]]
* [[Adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
* [[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 ==
* 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/#!NumericalMethods
* http://e6.ijs.si/ParallelAndDistributedSystems/#!utils
* http://e6.ijs.si/ParallelAndDistributedSystems/#!NUMA

Non-Newtonian fluid

2020-10-21T11:29:38Z

MihaR: Improved grammar

Click here to return back to [[Fluid Mechanics]]

=Introduction=
Other fluid mechanics examples have assumed that the fluids in question are Newtonian, but that is not always the case when dealing with real world applications. Many significant fluids exhibit non-Newtonian behaviour with examples including blood, foodstuffs, polymer melts... Newtonian fluids have a linear relationship between viscous stress and shear rate but that is no longer true for non-Newtonian fluids.

Some fluids exhibit increasing or decreasing viscosity with increasing shear rate which is the simplest form of non-Newtonian behaviour that we will tackle in the following sections. Some act as solids at small stress but begin flowing when sufficiently disturbed, while others behaviour depends on the duration of movement and applied stress.

=Model formulation=
Viscous forces are a result of velocity differentials within the fluid and are therefore connected with the velocity gradient. We are interested in the symmetric part that excludes the rotational contributions. Shear rate tensor $\dot{u}$ for incompressible fluids can be expressed as
$$
\dot{u}=\frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{v}+(\boldsymbol{\nabla} \boldsymbol{v})^{\top}\right),
$$
where $\boldsymbol{v}$ denotes the fluid velocity. We introduce viscosity $\mu$ through a generalized Newton's law of viscosity with viscous stress tensor $\tau$ expressed as
$$
\tau = \mu \dot{u}.
$$
In general $\mu$ is an arbitrary tensor but usually reduces to a scalar value due to isotropic properties of fluids. Non-Newtonian fluids are fluids where viscosity is a function of shear stress and duration of stress. There are many, often quantitative, models of varying complexities used to describe the viscosity of various non-Newtonian fluids.

We limit ourselves to Ostwald-de Waele power law model
$$
\mu = \mu_0 \left\| \dot{u} \right\| ^{n-1} = \mu_0 \left( 2 \dot{u}{:}\dot{u} \right) ^\frac{n-1}{2},
$$

<figure id="fig:devahl_scheme">
[[File:image.png|thumb|<caption>Scheme of the de Vahl Davis benchmark test.</caption>]]
</figure>

which is one of the simplest non-Newtonian viscosity models. We introduced two parameters used to fit the model to the experimental data for the fluid in question. Viscosity constant $\mu_0$ used as a scaling factor and non-Newtonian power index $n$ that describes the behaviour. Model with $n < 1$ describes shear-thinning (pseudoplastic), $n = 1$ Newtonian and $n > 1$ shear thickening (dilatant) fluids.

We will focus on shear thinning fluids as their behaviour of reduced viscosity at higher velocity gradients leads towards interesting behaviour with decreasing $n$.

=de Vahl Davis example=
We take a closer look at the behaviour of non-Newtonian fluid in a natural convection driven flow within a square cavity. The convection is driven by differentially heated east and west walls with the north and south walls thermally insulated as is schematically shown in <xr id="fig:devahl_scheme"/>.

This is a common benchmark with existing solutions that we can compare our meshless implementation against. Convergence and comparison with the benchmark article is shown in <xr id="fig:nonNewtonianConvergence"/> where we plot how the average Nusselt number changes with increasing node density. Nusselt number is a dimensionless number that expresses the ratio between the convective and conductive heat transfer. We use the average value on the western wall and define it as

<figure id="fig:nonNewtonianConvergence">
[[File:nonNewtonianConvergence.png|thumb|left|700px|
<caption>Convergence</caption> of the average Nusselt number value with the increasing node count and comparison with reference values from the reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref>. Convergence is calculated at Ra$=10^4$ and Pr$=10$.]]
</figure>

$$
\mathrm{Nu}=\frac{\Omega_W}{T_{H}-T_{C}}\left|\frac{\partial T}{\partial p_x}\right|_{p_x=0},
$$
where $T$ denotes the temperature, $p_x$ the $x$ coordinate and $\Omega_W$ the cavity width. Nusselt number is a convenient tool that can be used to reduce the dynamics of a complex multidimensional system into a single scalar value.

We define the problem in terms of Rayleigh and Prandtl dimensionless numbers that need to be slightly modified to remain relevant in a non-Newtonian case but their specifics are not really relevant for this review. Dimensionless number definitions used in these calculations match reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref> and can be viewed there.

Parameters that lead to an oscillatory behaviour in the de Vahl Davis cavity can be found in literature<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>. The referenced article found oscillatory behaviour in Newtonian fluids at Ra=$5 \cdot 10^4$ and Pr=$0.01$, but we want to use a slightly lower Ra value to find the transition as the dynamics become wilder with a smaller $n$. We chose Ra=$2 \cdot 10^4$ and Pr=$0.01$ as the parameters to be used for all subsequently displayed calculations.

The following figures display or mention dimensionless time that is defined with the same definition as described in the reference<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>.

<figure id="fig:nusseltEvolutionNN">
[[File:nusseltEvolutionNN.png|thumb|left|700px|
<caption>Evolution of the average Nusselt number in a differentially heated square cavity for different values of the non-Newtonian power index $n$.</caption>]]
</figure>

We reuse the Nusselt number to analyse the differences in behaviour for different non-Newtonian power indices $n$ shown in <xr id="fig:nusseltEvolutionNN"/>. The Newtonian case with $n=1$ settles into a stationary case which is expected, as we chose parameters below the previously identified oscillatory solution. The time evolutions of the average Nusselt number become oscillatory as we increase the strenght of non-Newtonian behaviour, while the solutions become less and less regular as we reduce the power index $n$, which can be seen in pictures of velocity and temperature distributions shown in <xr id="fig:squareCavityDisplay"/>. Detailed animation of behaviour in the most extreme $n=0.7$ case is shown in <xr id="fig:animation"/>.

Exact implementation of the non-Newtonian fluid modeled with Ostwald-de Waele power law can be found in the [[Non-Newtonian_fluid_example|example]].

<figure id="fig:squareCavityDisplay">
[[File:squareCavity_t10_Ra20000Pr001.png|frame|center|
<caption>Velocity and temperature fields within the cavity at $t=10$ for different values of the non-Newtonian power index $n$.</caption>]]
</figure>

<figure id="fig:animation">
[[File:squareCavityAnimationNN_n07.mp4|thumb|center|upright=2|
<caption>Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.</caption>]]
</figure>

=References=

<references />

Non-Newtonian fluid

2020-10-12T18:32:54Z

MihaR: Link to example

Click here to return back to [[Fluid Mechanics]]

=Introduction=
Other fluid mechanics examples have assumed that the fluids in question are Newtonian, but that is not always the case when dealing with real world applications. Many significant fluids exhibit non-Newtonian behaviour with examples including blood, foodstuffs, polymer melts... Newtonian fluids have a linear relationship between viscous stress and shear rate but that is no longer true for non-Newtonian fluids.

Some fluids exhibit increasing or decreasing viscosity with increasing shear rate which is the simplest form of non-Newtonian behaviour that we will tackle in the following sections. Some act as solids at small stress but begin flowing when sufficiently disturbed, while others behaviour depends on the duration of movement and applied stress.

=Model formulation=
Viscous forces are a result of velocity differentials within the fluid and are therefore connected with the velocity gradient. We are interested in the symmetric part that excludes the rotational contributions. Shear rate tensor $\dot{u}$ for incompressible fluids can be expressed as
$$
\dot{u}=\frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{v}+(\boldsymbol{\nabla} \boldsymbol{v})^{\top}\right),
$$
where $\boldsymbol{v}$ denotes the fluid velocity. We introduce viscosity $\mu$ through a generalized Newton's law of viscosity with viscous stress tensor $\tau$ expressed as
$$
\tau = \mu \dot{u}.
$$
In general $\mu$ is an arbitrary tensor but usually reduces to a scalar value due to isotropic properties of fluids. Non-Newtonian fluids are fluids where viscosity is a function of shear stress and duration of stress. There are many, often quantitative, models of varying complexities used to describe the viscosity of various non-Newtonian fluids.

We limit ourselves to Ostwald-de Waele power law model
$$
\mu = \mu_0 \left\| \dot{u} \right\| ^{n-1} = \mu_0 \left( 2 \dot{u}{:}\dot{u} \right) ^\frac{n-1}{2},
$$

<figure id="fig:devahl_scheme">
[[File:image.png|thumb|<caption>Scheme of the de Vahl Davis benchmark test.</caption>]]
</figure>

which is one of the simplest non-Newtonian viscosity models. We introduced two parameters used to fit the model to the experimental data for the fluid in question. Viscosity constant $\mu_0$ used as a scaling factor and non-Newtonian power index $n$ that describes the behaviour. Model with $n < 1$ describes shear-thinning (pseudoplastic), $n = 1$ Newtonian and $n > 1$ shear thickening (dilatant) fluids.

We will focus on shear thinning fluids as their behaviour of reduced viscosity at higher velocity gradients leads towards interesting behaviour with decreasing $n$.

=de Vahl Davis example=
We take a closer look at the behaviour of non-Newtonian fluid in a natural convection driven flow within a square cavity. The convection is driven by differentially heated east and west walls with the north and south walls thermally insulated as is schematically shown in <xr id="fig:devahl_scheme"/>.

This is a common benchmark with existing solutions that we can benchmark our meshless implementation against. Convergence and comparison with benchmark article is shown in <xr id="fig:nonNewtonianConvergence"/> where we plot the changes of average Nusselt number with increasing node density. Nusselt number is a dimensionless number that expresses the ratio between convective and conductive heat transfer. We use the average value on the western wall and define it as

<figure id="fig:nonNewtonianConvergence">
[[File:nonNewtonianConvergence.png|thumb|left|700px|
<caption>Convergence</caption> of average Nusselt number value with increasing node count and comparison with reference values from reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref>. Convergence is calculated at Ra$=10^4$ and Pr$=10$.]]
</figure>

$$
\mathrm{Nu}=\frac{\Omega_W}{T_{H}-T_{C}}\left|\frac{\partial T}{\partial p_x}\right|_{p_x=0},
$$
where $T$ denotest the demperature, $p_x$ $x$ coordinate and $\Omega_W$ cavity width. Nusselt number is a convenient tool that can be used to reduce the dynamics of a complex multidimensional system into a single scalar value.

We define the problem in terms of Rayleigh and Prandtl dimensionless numbers that need to be slightly modified to remain relevant in non-Newtonian case but specifics are not really relevant for this review. Dimensionless number definitions used in these calculations match reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref> and can be viewed there.

Parameters that lead to oscillatory behaviour in de Vahl Davis cavity can be found in literature<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>. Referenced article found oscillatory behaviour in Newtonian fluids at Ra=$5 \cdot 10^4$ and Pr=$0.01$ but we want to start slightly lower to find the transition as the dynamics become wilder with smaller $n$. We chose Ra=$2 \cdot 10^4$ and Pr=$0.01$ as the parameters to be used for all subsequently displayed calculations.

The following figures display or mention dimensionless time that is defined with the same definition as described in reference<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>.

<figure id="fig:nusseltEvolutionNN">
[[File:nusseltEvolutionNN.png|thumb|left|700px|
<caption>Evolution of average Nusselt number in a differentially heated square cavity for different values of non-Newtonian power index $n$.</caption>]]
</figure>

We reuse the Nusselt number to analyse the differences in behaviour for different non-Newtonian power indices $n$ shown in <xr id="fig:nusseltEvolutionNN"/>. The Newtonian case with $n=1$ settles into a stationary case which is expected, as we chose parameters below the previously identified oscillatory solution. Evolutions of Nusselt number become oscillatory as we increase the strenght of non-Newtonian behaviour while the solutions become less and less regular as we reduce the power index $n$ which can be seen in pictures of velocity and temperature distributions shown in <xr id="fig:squareCavityDisplay"/>. Detailed animation of behaviour in the most extreme $n=0.7$ case is shown in <xr id="fig:animation"/>.

Exact implementation of non-Newtonian fluid modeled with Ostwald-de Waele power law can be found in [[Non-Newtonian_fluid_example|example]].

<figure id="fig:squareCavityDisplay">
[[File:squareCavity_t10_Ra20000Pr001.png|frame|center|
<caption>Velocity and temperature fields within the cavity at $t=10$ for different values of non-Newtonian power index $n$.</caption>]]
</figure>

<figure id="fig:animation">
[[File:squareCavityAnimationNN_n07.mp4|thumb|center|upright=2|
<caption>Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.</caption>]]
</figure>

=References=

<references />

Non-Newtonian fluid

2020-10-12T18:13:54Z

MihaR: Added <caption> tags to figure captions

Click here to return back to [[Fluid Mechanics]]

=Introduction=
Other fluid mechanics examples have assumed that the fluids in question are Newtonian, but that is not always the case when dealing with real world applications. Many significant fluids exhibit non-Newtonian behaviour with examples including blood, foodstuffs, polymer melts... Newtonian fluids have a linear relationship between viscous stress and shear rate but that is no longer true for non-Newtonian fluids.

Some fluids exhibit increasing or decreasing viscosity with increasing shear rate which is the simplest form of non-Newtonian behaviour that we will tackle in the following sections. Some act as solids at small stress but begin flowing when sufficiently disturbed, while others behaviour depends on the duration of movement and applied stress.

=Model formulation=
Viscous forces are a result of velocity differentials within the fluid and are therefore connected with the velocity gradient. We are interested in the symmetric part that excludes the rotational contributions. Shear rate tensor $\dot{u}$ for incompressible fluids can be expressed as
$$
\dot{u}=\frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{v}+(\boldsymbol{\nabla} \boldsymbol{v})^{\top}\right),
$$
where $\boldsymbol{v}$ denotes the fluid velocity. We introduce viscosity $\mu$ through a generalized Newton's law of viscosity with viscous stress tensor $\tau$ expressed as
$$
\tau = \mu \dot{u}.
$$
In general $\mu$ is an arbitrary tensor but usually reduces to a scalar value due to isotropic properties of fluids. Non-Newtonian fluids are fluids where viscosity is a function of shear stress and duration of stress. There are many, often quantitative, models of varying complexities used to describe the viscosity of various non-Newtonian fluids.

We limit ourselves to Ostwald-de Waele power law model
$$
\mu = \mu_0 \left\| \dot{u} \right\| ^{n-1} = \mu_0 \left( 2 \dot{u}{:}\dot{u} \right) ^\frac{n-1}{2},
$$

<figure id="fig:devahl_scheme">
[[File:image.png|thumb|<caption>Scheme of the de Vahl Davis benchmark test.</caption>]]
</figure>

which is one of the simplest non-Newtonian viscosity models. We introduced two parameters used to fit the model to the experimental data for the fluid in question. Viscosity constant $\mu_0$ used as a scaling factor and non-Newtonian power index $n$ that describes the behaviour. Model with $n < 1$ describes shear-thinning (pseudoplastic), $n = 1$ Newtonian and $n > 1$ shear thickening (dilatant) fluids.

We will focus on shear thinning fluids as their behaviour of reduced viscosity at higher velocity gradients leads towards interesting behaviour with decreasing $n$.

=de Vahl Davis example=
We take a closer look at the behaviour of non-Newtonian fluid in a natural convection driven flow within a square cavity. The convection is driven by differentially heated east and west walls with the north and south walls thermally insulated as is schematically shown in <xr id="fig:devahl_scheme"/>.

This is a common benchmark with existing solutions that we can benchmark our meshless implementation against. Convergence and comparison with benchmark article is shown in <xr id="fig:nonNewtonianConvergence"/> where we plot the changes of average Nusselt number with increasing node density. Nusselt number is a dimensionless number that expresses the ratio between convective and conductive heat transfer. We use the average value on the western wall and define it as

<figure id="fig:nonNewtonianConvergence">
[[File:nonNewtonianConvergence.png|thumb|left|700px|
<caption>Convergence</caption> of average Nusselt number value with increasing node count and comparison with reference values from reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref>. Convergence is calculated at Ra$=10^4$ and Pr$=10$.]]
</figure>

$$
\mathrm{Nu}=\frac{\Omega_W}{T_{H}-T_{C}}\left|\frac{\partial T}{\partial p_x}\right|_{p_x=0},
$$
where $T$ denotest the demperature, $p_x$ $x$ coordinate and $\Omega_W$ cavity width. Nusselt number is a convenient tool that can be used to reduce the dynamics of a complex multidimensional system into a single scalar value.

We define the problem in terms of Rayleigh and Prandtl dimensionless numbers that need to be slightly modified to remain relevant in non-Newtonian case but specifics are not really relevant for this review. Dimensionless number definitions used in these calculations match reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref> and can be viewed there.

Parameters that lead to oscillatory behaviour in de Vahl Davis cavity can be found in literature<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>. Referenced article found oscillatory behaviour in Newtonian fluids at Ra=$5 \cdot 10^4$ and Pr=$0.01$ but we want to start slightly lower to find the transition as the dynamics become wilder with smaller $n$. We chose Ra=$2 \cdot 10^4$ and Pr=$0.01$ as the parameters to be used for all subsequently displayed calculations.

The following figures display or mention dimensionless time that is defined with the same definition as described in reference<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>.

<figure id="fig:nusseltEvolutionNN">
[[File:nusseltEvolutionNN.png|thumb|left|700px|
<caption>Evolution of average Nusselt number in a differentially heated square cavity for different values of non-Newtonian power index $n$.</caption>]]
</figure>

We reuse the Nusselt number to analyse the differences in behaviour for different non-Newtonian power indices $n$ shown in <xr id="fig:nusseltEvolutionNN"/>. The Newtonian case with $n=1$ settles into a stationary case which is expected, as we chose parameters below the previously identified oscillatory solution. Evolutions of Nusselt number become oscillatory as we increase the strenght of non-Newtonian behaviour while the solutions become less and less regular as we reduce the power index $n$ which can be seen in pictures of velocity and temperature distributions shown in <xr id="fig:squareCavityDisplay"/>. Detailed animation of behaviour in the most extreme $n=0.7$ case is shown in <xr id="fig:animation"/>.

Exact implementation of non-Newtonian fluid modeled with Ostwald-de Waele power law can be found in example.
TODO: link to example

<figure id="fig:squareCavityDisplay">
[[File:squareCavity_t10_Ra20000Pr001.png|frame|center|
<caption>Velocity and temperature fields within the cavity at $t=10$ for different values of non-Newtonian power index $n$.</caption>]]
</figure>

<figure id="fig:animation">
[[File:squareCavityAnimationNN_n07.mp4|thumb|center|upright=2|
<caption>Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.</caption>]]
</figure>

=References=

<references />

Cahn-Hilliard equation

2020-10-12T18:08:24Z

MihaR:

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

In this example we will show how to solve the [https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation Cahn-Hilliard] equation in Medusa. The solution provides an implementation of periodic boundary conditions and shows how to combine implicit and explicit operators. This example uses a custom biharmonic operator, so we recommend firstly reading the tutorial on [[Customization]].

== Cahn-Hilliard equation ==
The Cahn-Hilliard equation describes the separation of binary fluid into pure domains. We define concentration, a scalar field $c(x, y) \in [-1, 1]$ to describe the fluid, with values $c \pm 1$ representing pure phases. In this case the equation can be written as

\[
\frac{\partial c}{\partial t} = D \nabla^2 (c^3 - c - \gamma \nabla^2 c),
\]
where $D$ is the diffusion constant and $\sqrt{\gamma}$ is the characteristic transition length between domains.

We want to avoid any boundary effects and solve the equation for infinite fluid which can be approximated by using periodic boundary conditions.

== Periodic boundary conditions ==

<figure id="fig:periodicBands">
[[File:PeriodicBands.png|thumb|upright=2|<caption>Periodic band display.

Left subplot shows an example of filled domain with $border = 10$ and $dx = 0.2$. Interior nodes and their periodic counterparts share colors to illustrate periodicity.

Right subplot shows the neighbourhood of bottom left corner node $(-border, -border)$, with its support nodes highlighted in red. This example uses same fill parameters as the left figure and a support size of 100.
</caption>]]
</figure>

The simplest method of implementing a periodic boundary is adding a band of periodic nodes that mirror values from nodes that lie close to the opposite boundary. We have to create periodic nodes based on locations of interior nodes and construct a mapping that will allow us to update their values based on their interior counterparts. Width of the periodic band should be chosen with support size in mind. Periodic band should be wide enough to include all nodes required to calculate interior values. An example that can help us determine the required width of the periodic band is shown on the right subplot of <xr id="fig:periodicBands"/> and shows that we need a periodic band width of $6 dx$ for a support size of 100.

We start by creating a square domain centered on 0 and uniformly filling it with nodes. $border$ parameter determines the square size while $dx$, the typical distance between nodes, determines the node density.

<syntaxhighlight lang="cpp">
mm::BoxShape<vec_t> box(-border, border);
mm::DomainDiscretization<vec_t> domain = box.discretizeBoundaryWithStep(dx);
mm::GeneralFill<vec_t> fill;
fill.seed(randomSeed);
domain.fill(fill, dx);
</syntaxhighlight>

Afterwards we remove boundary nodes at positive $border$ and move others into interior. This creates a $[-border, border) \times [-border, border)$ interior that can be made periodic.

<syntaxhighlight lang="cpp">
Eigen::Matrix2d borders; // col1: - border, col2: + border, rows for dimensions
borders << domain.shape().bbox().first, domain.shape().bbox().second;
for (auto idx : domain.boundary()) {
const auto& pos = domain.pos(idx);
if (pos[0] != borders(0, 1) && pos[1] != borders(1, 1)) {
domain.changeToInterior(idx, pos, 1);
}
}
domain.removeBoundaryNodes();
</syntaxhighlight>

Now we are ready to find all nodes that lie sufficiently close to the border and create their periodic counterparts. In general rectangle case we would need to check all dimesions and create periodic nodes in all subsets of dimensions where that node lies close to the border. In this 2D case that simplifies to checking whether node lies close to the border in $x$, $y$ or both, which would mean a periodic corner. We could simplify our code and ignore the small periodic corners but that would lead to slightly worse results.

Periodic boundary nodes are created with type -2 so that we can still differentiate between them and any legitimate boundary nodes we would have in case of mixing periodic and normal boundary conditions. Boundary normal can be chosen arbitrarily, periodic nodes are only used as value holders.

During this process we also store periodic node indices in $periodicNodes$ and indices of the original nodes in $interiorMapping$.

<syntaxhighlight lang="cpp">
mm::Range<int> periodicNodes;
mm::Range<int> interiorMapping;
mm::Vec2d normal(0, 0); // not used in this case and can be chosen arbitrarily
for (auto idx : domain.interior()) {
const auto& pos = domain.pos(idx);
// col1: distance to - border, col2: distance to + border, rows represent dimensions
auto borderDistance = borders.colwise() - pos;
const auto& isPeriodic = borderDistance.array().abs() < periodicBandWidth * dx;
for (int i0 = 0; i0 < 2; i0++) {
// check for periodicity in 1st dimension (x)
if (isPeriodic(0, i0)) {
auto periodicPos = pos;
periodicPos(0) = borders(0, (i0 + 1) % 2) - borderDistance(0, i0);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
for (int i1 = 0; i1 < 2; i1++) {
// check for periodicity in 2nd dimension (y)
if (i0 == 0 && isPeriodic(1, i1)) {
auto periodicPos = pos;
periodicPos(1) = borders(1, (i1 + 1) % 2) - borderDistance(1, i1);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
// check for periodicity in both dimensions (corners)
if (isPeriodic(0, i0) && isPeriodic(1, i1)) {
auto periodicPos = pos;
periodicPos(0) = borders(0, (i0 + 1) % 2) - borderDistance(0, i0);
periodicPos(1) = borders(1, (i1 + 1) % 2) - borderDistance(1, i1);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
}
}
}
</syntaxhighlight>

== Solution ==
<figure id="fig:timeEvolution">
[[File:TimeEvolution.mp4|thumb|upright=2|<caption>Time evolution of the concentration. Initially mixed fluid separates into domains that grow with time.
</caption>]]
</figure>
We need to handle two uncommon terms to solve the equation. Biharmonic ($\nabla^4$) operator can be implemented with custom operators as shown in [[Customization]] example. Nonlinear ($\nabla^2 c^3$) term proves to be more problematic, as it can not be solved implicitly. We either need to solve the entire equation with explicit time stepping or include the nonlinear term in the right hand side of our system of equations.

Initial concentration is random.

=== Explicit solution ===
Explicit solution is simpler to implement and faster when we want to solve the equation for very small time steps.

We use Euler method for time stepping and copy interior node values to their periodic counterparts at every step.
<syntaxhighlight lang="cpp">
double t = 0;
Eigen::VectorXd bracket;
while (t <= tEnd) {
bracket = concentration.array() * concentration.array() * concentration.array() - concentration.array();
for (int i : allNodes) {
bracket[i] -= L * expOp.lap(concentration, i);
}
for (int i : interiorNodes) {
concentration[i] = concentration[i] + dt * D * expOp.lap(bracket, i);
}
t += dt;

for (int i = 0; i < periodicSize; i++) {
concentration[periodicNodes[i]] = concentration[interiorMapping[i]];
}
}
</syntaxhighlight>

<figure id="fig:timeEvolutionStatic">
[[File:cahnHilliardTimeEvolution.png|thumb|upright=1.5|<caption>Time evolution of the concentration. Initially mixed fluid separates into domains that grow with time.
</caption>]]
</figure>

=== Mixed solution ===
Mixed solution allows for much longer (several orders) time steps as we solve the majority of our equation implicitly. Explicit computation of the nonlinear term will always lead to inaccuracies, especially for long $dt$.

We set the equation on interior and identities on boundary.
<syntaxhighlight lang="cpp">
for (int i : interiorNodes) {
impOp.value(i) + D * dt * impOp.lap(i) + D * dt * L * impOp.apply<Biharmonic<dim>>(i) = rhs[i];
}
for (int i = 0; i < periodicSize; i++) {
impOp.value(periodicNodes[i]) = rhs[interiorMapping[i]];
}

Eigen::SparseLU<decltype(M)> solver;
solver.compute(M);
</syntaxhighlight>

We have to explicitly compute nonlinear $\nabla^2 c^3$ term at every step and add it to the rhs. Once we solve the system for new concentration we have to copy interior node values to their periodic counterparts.

<syntaxhighlight lang="cpp">
/// Time stepping.
double t = 0;
Eigen::VectorXd cubedConcentration;
while (t <= tEnd) {
/// Prepare rhs.
cubedConcentration = concentration.array() * concentration.array() * concentration.array();
for (int i : allNodes) {
concentration[i] = D * dt * expOp.lap(cubedConcentration, i) + concentration[i];
}

concentration = solver.solve(concentration);
t += dt;

/// Copy values from interior nodes to their periodic representations.
for (int i = 0; i < periodicSize; i++) {
concentration[periodicNodes[i]] = concentration[interiorMapping[i]];
}
}

</syntaxhighlight>

The complete example is avalible on [https://gitlab.com/e62Lab/medusa/blob/dev/examples/cahnHilliard_equation git].

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

File:CahnHilliardTimeEvolution.png

2020-10-12T18:01:40Z

MihaR: Example visualization of concentration evolution for Cahn-Hilliard equation.

Example visualization of concentration evolution for Cahn-Hilliard equation.

Non-Newtonian fluid

2020-10-08T19:42:27Z

MihaR: Animation display fix

Click here to return back to [[Fluid Mechanics]]

=Introduction=
Other fluid mechanics examples have assumed that the fluids in question are Newtonian, but that is not always the case when dealing with real world applications. Many significant fluids exhibit non-Newtonian behaviour with examples including blood, foodstuffs, polymer melts... Newtonian fluids have a linear relationship between viscous stress and shear rate but that is no longer true for non-Newtonian fluids.

Some fluids exhibit increasing or decreasing viscosity with increasing shear rate which is the simplest form of non-Newtonian behaviour that we will tackle in the following sections. Some act as solids at small stress but begin flowing when sufficiently disturbed, while others behaviour depends on the duration of movement and applied stress.

=Model formulation=
Viscous forces are a result of velocity differentials within the fluid and are therefore connected with the velocity gradient. We are interested in the symmetric part that excludes the rotational contributions. Shear rate tensor $\dot{u}$ for incompressible fluids can be expressed as
$$
\dot{u}=\frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{v}+(\boldsymbol{\nabla} \boldsymbol{v})^{\top}\right),
$$
where $\boldsymbol{v}$ denotes the fluid velocity. We introduce viscosity $\mu$ through a generalized Newton's law of viscosity with viscous stress tensor $\tau$ expressed as
$$
\tau = \mu \dot{u}.
$$
In general $\mu$ is an arbitrary tensor but usually reduces to a scalar value due to isotropic properties of fluids. Non-Newtonian fluids are fluids where viscosity is a function of shear stress and duration of stress. There are many, often quantitative, models of varying complexities used to describe the viscosity of various non-Newtonian fluids.

We limit ourselves to Ostwald-de Waele power law model
$$
\mu = \mu_0 \left\| \dot{u} \right\| ^{n-1} = \mu_0 \left( 2 \dot{u}{:}\dot{u} \right) ^\frac{n-1}{2},
$$

<figure id="fig:devahl_scheme">
[[File:image.png|thumb|Scheme of the de Vahl Davis benchmark test.]]
</figure>

which is one of the simplest non-Newtonian viscosity models. We introduced two parameters used to fit the model to the experimental data for the fluid in question. Viscosity constant $\mu_0$ used as a scaling factor and non-Newtonian power index $n$ that describes the behaviour. Model with $n < 1$ describes shear-thinning (pseudoplastic), $n = 1$ Newtonian and $n > 1$ shear thickening (dilatant) fluids.

We will focus on shear thinning fluids as their behaviour of reduced viscosity at higher velocity gradients leads towards interesting behaviour with decreasing $n$.

=de Vahl Davis example=
We take a closer look at the behaviour of non-Newtonian fluid in a natural convection driven flow within a square cavity. The convection is driven by differentially heated east and west walls with the north and south walls thermally insulated as is schematically shown in <xr id="fig:devahl_scheme"/>.

This is a common benchmark with existing solutions that we can benchmark our meshless implementation against. Convergence and comparison with benchmark article is shown in <xr id="fig:nonNewtonianConvergence"/> where we plot the changes of average Nusselt number with increasing node density. Nusselt number is a dimensionless number that expresses the ratio between convective and conductive heat transfer. We use the average value on the western wall and define it as

<figure id="fig:nonNewtonianConvergence">
[[File:nonNewtonianConvergence.png|thumb|left|700px|
Convergence of average Nusselt number value with increasing node count and comparison with reference values from reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref>. Convergence is calculated at Ra$=10^4$ and Pr$=10$.]]
</figure>

$$
\mathrm{Nu}=\frac{\Omega_W}{T_{H}-T_{C}}\left|\frac{\partial T}{\partial p_x}\right|_{p_x=0},
$$
where $T$ denotest the demperature, $p_x$ $x$ coordinate and $\Omega_W$ cavity width. Nusselt number is a convenient tool that can be used to reduce the dynamics of a complex multidimensional system into a single scalar value.

We define the problem in terms of Rayleigh and Prandtl dimensionless numbers that need to be slightly modified to remain relevant in non-Newtonian case but specifics are not really relevant for this review. Dimensionless number definitions used in these calculations match reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref> and can be viewed there.

Parameters that lead to oscillatory behaviour in de Vahl Davis cavity can be found in literature<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>. Referenced article found oscillatory behaviour in Newtonian fluids at Ra=$5 \cdot 10^4$ and Pr=$0.01$ but we want to start slightly lower to find the transition as the dynamics become wilder with smaller $n$. We chose Ra=$2 \cdot 10^4$ and Pr=$0.01$ as the parameters to be used for all subsequently displayed calculations.

The following figures display or mention dimensionless time that is defined with the same definition as described in reference<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>.

<figure id="fig:nusseltEvolutionNN">
[[File:nusseltEvolutionNN.png|thumb|left|700px|
Evolution of average Nusselt number in a differentially heated square cavity for different values of non-Newtonian power index $n$.]]
</figure>

We reuse the Nusselt number to analyse the differences in behaviour for different non-Newtonian power indices $n$ shown in <xr id="fig:nusseltEvolutionNN"/>. The Newtonian case with $n=1$ settles into a stationary case which is expected, as we chose parameters below the previously identified oscillatory solution. Evolutions of Nusselt number become oscillatory as we increase the strenght of non-Newtonian behaviour while the solutions become less and less regular as we reduce the power index $n$ which can be seen in pictures of velocity and temperature distributions shown in <xr id="fig:squareCavityDisplay"/>. Detailed animation of behaviour in the most extreme $n=0.7$ case is shown in <xr id="fig:animation"/>.

Exact implementation of non-Newtonian fluid modeled with Ostwald-de Waele power law can be found in example.
TODO: link to example

<figure id="fig:squareCavityDisplay">
[[File:squareCavity_t10_Ra20000Pr001.png|frame|center|
Velocity and temperature fields within the cavity at $t=10$ for different values of non-Newtonian power index $n$.]]
</figure>

<figure id="fig:animation">
[[File:squareCavityAnimationNN_n07.mp4|thumb|center|upright=2|
Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.]]
</figure>

=References=

<references />

Non-Newtonian fluid

2020-10-08T19:37:35Z

MihaR: Initial version

Click here to return back to [[Fluid Mechanics]]

=Introduction=
Other fluid mechanics examples have assumed that the fluids in question are Newtonian, but that is not always the case when dealing with real world applications. Many significant fluids exhibit non-Newtonian behaviour with examples including blood, foodstuffs, polymer melts... Newtonian fluids have a linear relationship between viscous stress and shear rate but that is no longer true for non-Newtonian fluids.

Some fluids exhibit increasing or decreasing viscosity with increasing shear rate which is the simplest form of non-Newtonian behaviour that we will tackle in the following sections. Some act as solids at small stress but begin flowing when sufficiently disturbed, while others behaviour depends on the duration of movement and applied stress.

=Model formulation=
Viscous forces are a result of velocity differentials within the fluid and are therefore connected with the velocity gradient. We are interested in the symmetric part that excludes the rotational contributions. Shear rate tensor $\dot{u}$ for incompressible fluids can be expressed as
$$
\dot{u}=\frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{v}+(\boldsymbol{\nabla} \boldsymbol{v})^{\top}\right),
$$
where $\boldsymbol{v}$ denotes the fluid velocity. We introduce viscosity $\mu$ through a generalized Newton's law of viscosity with viscous stress tensor $\tau$ expressed as
$$
\tau = \mu \dot{u}.
$$
In general $\mu$ is an arbitrary tensor but usually reduces to a scalar value due to isotropic properties of fluids. Non-Newtonian fluids are fluids where viscosity is a function of shear stress and duration of stress. There are many, often quantitative, models of varying complexities used to describe the viscosity of various non-Newtonian fluids.

We limit ourselves to Ostwald-de Waele power law model
$$
\mu = \mu_0 \left\| \dot{u} \right\| ^{n-1} = \mu_0 \left( 2 \dot{u}{:}\dot{u} \right) ^\frac{n-1}{2},
$$

<figure id="fig:devahl_scheme">
[[File:image.png|thumb|Scheme of the de Vahl Davis benchmark test.]]
</figure>

which is one of the simplest non-Newtonian viscosity models. We introduced two parameters used to fit the model to the experimental data for the fluid in question. Viscosity constant $\mu_0$ used as a scaling factor and non-Newtonian power index $n$ that describes the behaviour. Model with $n < 1$ describes shear-thinning (pseudoplastic), $n = 1$ Newtonian and $n > 1$ shear thickening (dilatant) fluids.

We will focus on shear thinning fluids as their behaviour of reduced viscosity at higher velocity gradients leads towards interesting behaviour with decreasing $n$.

=de Vahl Davis example=
We take a closer look at the behaviour of non-Newtonian fluid in a natural convection driven flow within a square cavity. The convection is driven by differentially heated east and west walls with the north and south walls thermally insulated as is schematically shown in <xr id="fig:devahl_scheme"/>.

This is a common benchmark with existing solutions that we can benchmark our meshless implementation against. Convergence and comparison with benchmark article is shown in <xr id="fig:nonNewtonianConvergence"/> where we plot the changes of average Nusselt number with increasing node density. Nusselt number is a dimensionless number that expresses the ratio between convective and conductive heat transfer. We use the average value on the western wall and define it as

<figure id="fig:nonNewtonianConvergence">
[[File:nonNewtonianConvergence.png|thumb|left|700px|
Convergence of average Nusselt number value with increasing node count and comparison with reference values from reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref>. Convergence is calculated at Ra$=10^4$ and Pr$=10$.]]
</figure>

$$
\mathrm{Nu}=\frac{\Omega_W}{T_{H}-T_{C}}\left|\frac{\partial T}{\partial p_x}\right|_{p_x=0},
$$
where $T$ denotest the demperature, $p_x$ $x$ coordinate and $\Omega_W$ cavity width. Nusselt number is a convenient tool that can be used to reduce the dynamics of a complex multidimensional system into a single scalar value.

We define the problem in terms of Rayleigh and Prandtl dimensionless numbers that need to be slightly modified to remain relevant in non-Newtonian case but specifics are not really relevant for this review. Dimensionless number definitions used in these calculations match reference<ref name="Turan2011">O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law Fluids in a square enclosure with differentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).</ref> and can be viewed there.

Parameters that lead to oscillatory behaviour in de Vahl Davis cavity can be found in literature<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>. Referenced article found oscillatory behaviour in Newtonian fluids at Ra=$5 \cdot 10^4$ and Pr=$0.01$ but we want to start slightly lower to find the transition as the dynamics become wilder with smaller $n$. We chose Ra=$2 \cdot 10^4$ and Pr=$0.01$ as the parameters to be used for all subsequently displayed calculations.

The following figures display or mention dimensionless time that is defined with the same definition as described in reference<ref name="Kosec2011">G. Kosec in B. Šarler, Solution of a low Prandtl number natural convection benchmark by a local meshless method, International Journal of Numerical Methods for Heat & Fluid Flow 23, 22 (2013).</ref>.

<figure id="fig:nusseltEvolutionNN">
[[File:nusseltEvolutionNN.png|thumb|left|700px|
Evolution of average Nusselt number in a differentially heated square cavity for different values of non-Newtonian power index $n$.]]
</figure>

We reuse the Nusselt number to analyse the differences in behaviour for different non-Newtonian power indices $n$ shown in <xr id="fig:nusseltEvolutionNN"/>. The Newtonian case with $n=1$ settles into a stationary case which is expected, as we chose parameters below the previously identified oscillatory solution. Evolutions of Nusselt number become oscillatory as we increase the strenght of non-Newtonian behaviour while the solutions become less and less regular as we reduce the power index $n$ which can be seen in pictures of velocity and temperature distributions shown in <xr id="fig:squareCavityDisplay"/>. Detailed animation of behaviour in the most extreme $n=0.7$ case is shown in <xr id="fig:animation"/>.

Exact implementation of non-Newtonian fluid modeled with Ostwald-de Waele power law can be found in example.
TODO: link to example

<figure id="fig:squareCavityDisplay">
[[File:squareCavity_t10_Ra20000Pr001.png|frame|center|
Velocity and temperature fields within the cavity at $t=10$ for different values of non-Newtonian power index $n$.]]
</figure>

<figure id="fig:animation">
[[File:squareCavityAnimationNN_n07.mp4|frame|center|
Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.]]
</figure>

=References=

<references />

File:SquareCavityAnimationNN n07.mp4

2020-10-08T19:34:35Z

MihaR: Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.

Animation of velocity and temperature fields within the cavity between $t=10$ and $t=11$ for $n=0.7$.

File:SquareCavity t10 Ra20000Pr001.png

2020-10-08T18:43:24Z

MihaR: Velocity magnitude and temperature within the cavity for different values of non-Newtonian power index $n$.

Velocity magnitude and temperature within the cavity for different values of non-Newtonian power index $n$.

File:NusseltEvolutionNN.png

2020-10-08T18:41:20Z

MihaR: Evolution of average Nusselt number in a differentially heated square cavity for different values of non-Newtonian power index $n$.

Evolution of average Nusselt number in a differentially heated square cavity for different values of non-Newtonian power index $n$.

File:NonNewtonianConvergence.png

2020-10-08T17:39:46Z

MihaR: MihaR uploaded a new version of File:NonNewtonianConvergence.png

Convergence of average Nusselt number value with increasing node count and comparison with reference values from O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law �uids in a square enclosure with di�erentially heated side walls subjected to constant temperatures, Journal of Non-Newtonian Fluid Mechanics 166, 1049 (2011).

File:NonNewtonianConvergence.png

2020-10-08T17:03:57Z

MihaR: Convergence of average Nusselt number value with increasing node count and comparison with reference values from O. Turan, A. Sachdeva, N. Chakraborty in R. J. Poole, Laminar natural convection of power-law �uids in a square enclosure with di�erent...

File:NonNewtonianConvergence.pdf

2020-10-08T16:53:03Z

MihaR:

De Vahl Davis natural convection test

2020-10-07T11:16:25Z

MihaR: Dead link fix

Click here to return back to [[Fluid Mechanics]]

=Intro=
The classical de Vahl Davis benchmark test is defined for the natural convection of the air ($\Pr =0.71$) in the square closed cavity (${{\text{A}}_{\text{R}}}=1$). The only physical free parameter of the test remains the thermal Rayleigh number. In the original paper [1] de Vahl Davis tested the problem up to the Rayleigh number ${{10}^{6}}$, however in the latter publications, the results of more intense simulations were presented with the Rayleigh number up to ${{10}^{8}}$. Lage and Bejan [2] showed that the laminar domain of the closed cavity natural convection problem is roughly below $\text{Gr1}{{\text{0}}^{9}}$. It was reported [3, 4] that the natural convection becomes unsteady for $\text{Ra}=2\cdot {{10}^5}$. Here we present a MLSM solution of the case.

\begin{equation}
\text{Ra}\text{=}\,\frac{\left| \mathbf{g} \right|{{\beta }_{T}}\left( {{T}_{H}}-{{T}_{C}} \right){{\Omega }_{H}}^{3}{{\rho }^{2}}{{c}_{p}}}{\lambda \mu }
\end{equation}
\begin{equation}
\text{Pr}=\frac{\mu {{c}_{p}}}{\lambda }
\end{equation}

[1] de Vahl Davis G. Natural convection of air in a square cavity: a bench mark numerical solution. Int J Numer Meth Fl. 1983;3:249-64.

[2] Lage JL, Bejan A. The Ra-Pr domain of laminar natural convection in an enclosure heated from the side. Numer Heat Transfer. 1991;A19:21-41.

[3] Janssen RJA, Henkes RAWM. Accuracy of finite-volume disretizations for the bifurcating natural-convection flow in a square cavity. Numer Heat Transfer. 1993;B24:191-207.

[4] Nobile E. Simulation of time-dependent flow in cavities with the additive-correction multigrid method, part II: Apllications. Numer Heat Transfer. 1996;B30:341-50.

<figure id="fig:devahl_scheme">
[[File:image.png]].
<caption>Scheme of the de Vahl Davis benchmark test </caption>
</figure>

=Code=
Full examples can be found under the examples in the [https://gitlab.com/e62Lab/medusa code repository].

==Explicit ACM method with CBS looks==
<syntaxhighlight lang="c++" line>
v2[boundary] = vec_t{0.0, 0.0};
T2[left] = O.T_cold;
T2[right] = O.T_hot;
//Time stepping
for (int step = 0; step <= O.t_steps; ++step) {
for (int i_count = 1; i_count < _MAX_ITER_; ++i_count) {
// Navier Stokes
for (auto c : interior) {
v2[c] = v1[c] + O.dt * (-1 / O.rho * op.grad(P1, c)
+ O.mu / O.rho * op.lap(v1, c)
- op.grad(v1, c) * v1[c]
+ O.g * (1 - O.beta * (T1[c] - O.T_ref)));
}

//Speed of sound
Range<scal_t> norm = v2.map([](const vec_t& p) { return p.norm(); });
scal_t C = O.dl * std::max(*std::max_element(norm.begin(), norm.end()), O.v_ref);
// Mass continuity
Range<scal_t> div_v;
for (auto c:all) {
div_v[c] = op.div(v2, c);
P2[c] = P1[c] - C * C * O.dt * O.rho * div_v[c] +
O.dl2 * C * C * O.dt * O.dt * op.lap(P1, c);
}
P1.swap(P2);
}

//heat transport
for (auto c : interior) {
T2[c] = T1[c] + O.dt * O.lam / O.rho / O.c_p * op.lap(T1, c) -
O.dt * v1[c].transpose() * op.grad(T1, c);
}
for (auto c : top) T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
for (auto c : bottom) T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
}
</syntaxhighlight>

== Explicit pressure correction ==
The solution of heat equation is the same as in above example
<syntaxhighlight lang="c++" line>
for (int step = 0; step <= O.t_steps; ++step) {
// Explicit Navier-Stokes computed on whole domain, including boundaries
// without pressure
for (int c:all) {
v_2[c] = v_1[c] + O.dt (
O.mu / O.rho * op.lap(v_1, c)
- op.grad(v_1, c) * v_1[c]
+ O.g * (1 - O.beta * (T_1[c] - O.T_ref)));
}
// Pressure correction
VecXd rhs_pressure(N + 1, 0); //Note N+1, +1 stands for regularization equation
rhs_pressure(N) = 0; // = 0 part of the regularization equation
for (int i:interior) rhs_pressure(c) = O.rho / O.dt * op.div(v_2, c);
for (int i: boundary) rhs_pressure(c) = O.rho / O.dt * v_2[c].dot(domain.normal(c));
VecXd solution = solver_p.solve(rhs_pressure);
alpha = solution[N];
VecXd P_c = solution.head(N);
for ( int i = interior) v_2[c] -= O.dt / O.rho * op.grad(P_c, c);
v_2[boundary] = 0; // force boundary conditions
//heat transport
for (auto c : interior) {
T2[c] = T1[c] + O.dt * O.lam / O.rho / O.c_p * op.lap(T1, c) -
O.dt * v1[c].transpose() * op.grad(T1, c);
}
for (auto c : top) T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
for (auto c : bottom) T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
}
</syntaxhighlight>

== Implicit ==
<syntaxhighlight lang="c++" line>
for (int step = 0; step <= O.t_steps_i; ++step) {
time_1 = std::chrono::high_resolution_clock::now();
// NAVIER STOKES
M_velocity = mat_t(2 * N, 2 * N);
// system
M_velocity.reserve(Range<int>(2 * N, O.n));
for (int i : interior) {
op.valuevec(M_velocity, i, 1 / O.dt_i);
op.lapvec(M_velocity, i, -O.mu / O.rho);
op.gradvec(M_velocity, i, v_1[i]);
}
for (int i : boundary) op.valuevec(M_velocity, i, 1); //sets velocity to 0

M_velocity.makeCompressed();
solver_v.compute(M_velocity);
// solution
Range <vec_t> rhs_vec(domain.size(), 0);
for (int i : interior) {
rhs_vec[i] = -1 / O.rho * op.grad(P, i) +
v_1[i] / O.dt_i
+ O.g * (1 - O.beta * (T[i] - O.T_ref));
}
// for (int i:top) rhs_vec[i] = vec_t{0,1};

v_2 = reshape<2>(solver_v.solveWithGuess(reshape(rhs_vec), reshape(v_1)));
// END OF NAVIER STOKES

// PRESSURE CORRECTION
VecXd rhs_pressure(N + 1, 0); //Note N+1, +1 stands for regularization equation
rhs_pressure(N) = 0; // = 0 part of the regularization equation
double alpha;
for (int i : interior) rhs_pressure(i) = O.rho / O.dt_i * op.div(v_2, i);
for (int i : boundary) rhs_pressure(i) = O.rho / O.dt * v_2[i].dot(domain.normal(i));

VecXd solution = solver_p.solve(rhs_pressure);
alpha = solution[N];
VecXd P_c = solution.head(N);
// apply velocity correction
for (int i : interior) {
v_2[i] -= O.dl * O.dt_i / O.rho * op.grad(P_c, i);
}
P += O.dl * P_c;
// enforce velocity BC
// v_2[boundary] = 0;
// END OF PRESSURE CORRECTION

// HEAT TRANSPORT
M_temperature = mat_t(N, N);
Range<int> per_row(N, O.n);
M_temperature.reserve(per_row);
// outer boundary dirichlet BC
for (int i : top) op.neumann_implicit(M_temperature, i, domain.normal(i), 1);
for (int i : bottom) op.neumann_implicit(M_temperature, i, domain.normal(i), 1);
for (int i : left) op.value(M_temperature, i, 1.0);
for (int i : right) op.value(M_temperature, i, 1.0);
// heat transport in air
for (int i : interior) {
op.value(M_temperature, i, 1.0); // time dependency
op.lap(M_temperature, i, -O.dt_i * O.lam / O.rho / O.c_p); //laplace in interior
op.grad(M_temperature, i, O.dt * v_2[i]);
}
M_temperature.makeCompressed();
solver_T.compute(M_temperature);

VectorXd rhs = VectorXd::Zero(N);
for (int i : interior) rhs(i) = T(i);
for (int i : top) rhs(i) = 0;
for (int i : bottom) rhs(i) = 0;
for (int i : left) rhs(i) = O.T_hot;
for (int i : right) rhs(i) = O.T_cold;
T = solver_T.solveWithGuess(rhs, T);
// END OF HEAT TRANSPORT
v_1.swap(v_2);
}
</syntaxhighlight>

=Results=
== Comparison of MLSM solution with reference data ==
Following video shows evolution of temperature and velocity magnitude for the $Ra=10^8$ case.

[[File:DeVahlDavisRa1e8.mp4|640px]]

In below galley you can find temperature contour plots, velocity magnitude contour plots, v_max and average hot side Nusselt number convergence behavior. The reference values are from:
*[a] de Vahl Davis G. Natural convection of air in a square cavity: a bench mark numerical solution. Int J Numer Meth Fl. 1983;3:249-64.
*[b] Sadat H, Couturier S. Performance and accuracy of a meshless method for laminar natural convection. Numer Heat Transfer. 2000;B37:455-67.
*[c] Wan DC, Patnaik BSV, Wei GW. A new benchmark quality solution for the buoyancy-driven cavity by discrete singular convolution. Numer Heat Transfer. 2001;B40:199-228.
*[d] Šarler B. A radial basis function collocation approach in computational fluid dynamics. CMES-Comp Model Eng. 2005;7:185-93.
*[e] Kosec G, Šarler B. Solution of thermo-fluid problems by collocation with local pressure correction. Int J Numer Method H. 2008;18:868-82.

<gallery>
File:DVD_cont_T_Ra3.png
File:DVD_cont_T_Ra4.png
File:DVD_cont_T_Ra5.png
File:DVD_cont_T_Ra6.png
File:DVD_cont_T_Ra7.png
File:DVD_cont_T_Ra8.png
File:DVD_cont_v_Ra3.png
File:DVD_cont_v_Ra4.png
File:DVD_cont_v_Ra5.png
File:DVD_cont_v_Ra6.png
File:DVD_cont_v_Ra7.png
File:DVD_cont_v_Ra8.png
File:DVD_Nu_convRa3.png
File:DVD_Nu_convRa4.png
File:DVD_Nu_convRa5.png
File:DVD_Nu_convRa6.png
File:DVD_Nu_convRa7.png
File:DVD_Nu_convRa8.png
File:DVD_vMax_convRa3.png
File:DVD_vMax_convRa4.png
File:DVD_vMax_convRa5.png
File:DVD_vMax_convRa6.png
File:DVD_vMax_convRa7.png
File:DVD_vMax_convRa8.png
</gallery>

Medusa

2020-03-23T20:40:33Z

MihaR: /* 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, 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]]
* [[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]]
* [[Customization | Customization example: Biharmonic equation]]
* [[Ghost nodes]]
* [[Electromagnetic scattering]]
* [[Schrödinger equation]]
* [[Wave equation]]
* [[Cahn-Hilliard equation]]
* [[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]]
* [[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]]

== 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.
* Basic PDE solutions
** [[Convection Diffusion equation | Convection Diffusion equation]]
**[[Wave equation application]]
* [[Adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
* [[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]]
* 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 ==
* 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/#!NumericalMethods
* http://e6.ijs.si/ParallelAndDistributedSystems/#!utils
* http://e6.ijs.si/ParallelAndDistributedSystems/#!NUMA

Bioheat equation

2020-03-23T20:27:37Z

MihaR:

<figure id="fig:brainBioheat">
[[File:BrainBioheat.png|thumb|upright=2|<caption>
</caption>]]
</figure>

The Pennes' bioheat equation is a standard model for temperature distrubution in living tissues, that enhances diffusion equation with a linear term describing blood flow and a constant metabolic heat source.

\[
\rho c \frac{\partial T}{\partial t} = \nabla(\lambda \nabla T) + W_b (T_a -T) + Q_m
\]

This example implements the stationary form of bioheat equation

\[
\nabla(\lambda \nabla T) + W_b (T_a -T) + Q_m = 0
\]
with Robin boundary conditions
\[
\lambda \frac{\partial T}{\partial \hat{n}} = h_s(T - T_{ext})
\]
on a human brain model
<ref name="Cvetkovic2016">
Mario Cvetković, Dragan Poljak, Akimasa Hirata, The electromagnetic-thermal dosimetry for the homogeneous human brain model, Engineering Analysis with Boundary Elements, Volume 63, 2016, Pages 61-73, ISSN 0955 7997, https://doi.org/10.1016/j.enganabound.2015.11.002.</ref>. Simulation nodes are based on FEM elements used in the referenced article with constants set to the default values from table 2 of the article.

Obtained solution is displayed on <xr id="fig:brainBioheat"/>

=References=
<references/>

Bioheat equation

2020-03-23T20:25:33Z

MihaR:

<figure id="fig:brainBioheat">
[[File:BrainBioheat.png|thumb|upright=2|<caption>
</caption>]]
</figure>

The Pennes' bioheat equation is a standard model for temperature distrubution in living tissues that enhances diffusion equation with a linear term describing blood flow and constant metabolic heat sources.

\[
\rho c \frac{\partial T}{\partial t} = \nabla(\lambda \nabla T) + W_b (T_a -T) + Q_m
\]

This example implements the stationary form of bioheat equation

\[
\nabla(\lambda \nabla T) + W_b (T_a -T) + Q_m = 0
\]
with Robin boundary conditions
\[
\lambda \frac{\partial T}{\partial \hat{n}} = h_s(T - T_{ext})
\]
on a human brain model
<ref name="Cvetkovic2016">
Mario Cvetković, Dragan Poljak, Akimasa Hirata, The electromagnetic-thermal dosimetry for the homogeneous human brain model, Engineering Analysis with Boundary Elements, Volume 63, 2016, Pages 61-73, ISSN 0955 7997, https://doi.org/10.1016/j.enganabound.2015.11.002.</ref>. Simulation nodes are based on the FEM elements used in the referenced article with constants set to the default values from table 2 of the article.

Obtained solution is displayed on <xr id="fig:brainBioheat"/>

=References=
<references/>

File:BrainBioheat.png

2020-03-23T20:24:51Z

MihaR:

Bioheat equation

2020-03-23T20:24:25Z

MihaR: Created page with "The Pennes' bioheat equation is a standard model for temperature distrubution in living tissues that enhances diffusion equation with a linear term describing blood flow and c..."

The Pennes' bioheat equation is a standard model for temperature distrubution in living tissues that enhances diffusion equation with a linear term describing blood flow and constant metabolic heat sources.

\[
\rho c \frac{\partial T}{\partial t} = \nabla(\lambda \nabla T) + W_b (T_a -T) + Q_m
\]

This example implements the stationary form of bioheat equation

\[
\nabla(\lambda \nabla T) + W_b (T_a -T) + Q_m = 0
\]
with Robin boundary conditions
\[
\lambda \frac{\partial T}{\partial \hat{n}} = h_s(T - T_{ext})
\]
on a human brain model
<ref name="Cvetkovic2016">
Mario Cvetković, Dragan Poljak, Akimasa Hirata, The electromagnetic-thermal dosimetry for the homogeneous human brain model, Engineering Analysis with Boundary Elements, Volume 63, 2016, Pages 61-73, ISSN 0955 7997, https://doi.org/10.1016/j.enganabound.2015.11.002.</ref>. Simulation nodes are based on the FEM elements used in the referenced article with constants set to the default values from table 2 of the article.

Obtained solution is displayed on <xr id="fig:brainBioheat"/>

<figure id="fig:brainBioheat">
[[File:BrainBioheat.png|thumb|upright=2|<caption>
</caption>]]
</figure>

=References=
<references/>

File:TimeEvolution.mp4

2020-02-13T11:40:39Z

MihaR: MihaR uploaded a new version of File:TimeEvolution.mp4

Medusa

2020-02-13T11:14:45Z

MihaR: /* Examples */

'''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, 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]]
* [[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]]
* [[Customization | Customization example: Biharmonic equation]]
* [[Ghost nodes]]
* [[Electromagnetic scattering]]
* [[Schrödinger equation]]
* [[Wave equation]]
* [[Cahn-Hilliard equation]]
* [[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]]
* [[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]]

== 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.
* Basic PDE solutions
** [[Convection Diffusion equation | Convection Diffusion equation]]
**[[Wave equation application]]
* [[Adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
* [[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]]
* 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>

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

Cahn-Hilliard equation

2020-02-11T14:10:53Z

MihaR:

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

In this example we will show how to solve the [https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation Cahn-Hilliard] equation in Medusa. The solution provides an implementation of periodic boundary conditions and shows how to combine implicit and explicit operators. This example uses a custom biharmonic operator, so we recommend firstly reading the tutorial on [[Customization]].

== Cahn-Hilliard equation ==
The Cahn-Hilliard equation describes the separation of binary fluid into pure domains. We define concentration, a scalar field $c(x, y) \in [-1, 1]$ to describe the fluid, with values $c \pm 1$ representing pure phases. In this case the equation can be written as

\[
\frac{\partial c}{\partial t} = D \nabla^2 (c^3 - c - \gamma \nabla^2 c),
\]
where $D$ is the diffusion constant and $\sqrt{\gamma}$ is the characteristic transition length between domains.

We want to avoid any boundary effects and solve the equation for infinite fluid which can be approximated by using periodic boundary conditions.

== Periodic boundary conditions ==

<figure id="fig:periodicBands">
[[File:PeriodicBands.png|thumb|upright=2|<caption>Periodic band display.

Left subplot shows an example of filled domain with $border = 10$ and $dx = 0.2$. Interior nodes and their periodic counterparts share colors to illustrate periodicity.

Right subplot shows the neighbourhood of bottom left corner node $(-border, -border)$, with its support nodes highlighted in red. This example uses same fill parameters as the left figure and a support size of 100.
</caption>]]
</figure>

The simplest method of implementing a periodic boundary is adding a band of periodic nodes that mirror values from nodes that lie close to the opposite boundary. We have to create periodic nodes based on locations of interior nodes and construct a mapping that will allow us to update their values based on their interior counterparts. Width of the periodic band should be chosen with support size in mind. Periodic band should be wide enough to include all nodes required to calculate interior values. An example that can help us determine the required width of the periodic band is shown on the right subplot of <xr id="fig:periodicBands"/> and shows that we need a periodic band width of $6 dx$ for a support size of 100.

We start by creating a square domain centered on 0 and uniformly filling it with nodes. $border$ parameter determines the square size while $dx$, the typical distance between nodes, determines the node density.

<syntaxhighlight lang="cpp">
mm::BoxShape<vec_t> box(-border, border);
mm::DomainDiscretization<vec_t> domain = box.discretizeBoundaryWithStep(dx);
mm::GeneralFill<vec_t> fill;
fill.seed(randomSeed);
domain.fill(fill, dx);
</syntaxhighlight>

Afterwards we remove boundary nodes at positive $border$ and move others into interior. This creates a $[-border, border) \times [-border, border)$ interior that can be made periodic.

<syntaxhighlight lang="cpp">
Eigen::Matrix2d borders; // col1: - border, col2: + border, rows for dimensions
borders << domain.shape().bbox().first, domain.shape().bbox().second;
for (auto idx : domain.boundary()) {
const auto& pos = domain.pos(idx);
if (pos[0] != borders(0, 1) && pos[1] != borders(1, 1)) {
domain.changeToInterior(idx, pos, 1);
}
}
domain.removeBoundaryNodes();
</syntaxhighlight>

Now we are ready to find all nodes that lie sufficiently close to the border and create their periodic counterparts. In general rectangle case we would need to check all dimesions and create periodic nodes in all subsets of dimensions where that node lies close to the border. In this 2D case that simplifies to checking whether node lies close to the border in $x$, $y$ or both, which would mean a periodic corner. We could simplify our code and ignore the small periodic corners but that would lead to slightly worse results.

Periodic boundary nodes are created with type -2 so that we can still differentiate between them and any legitimate boundary nodes we would have in case of mixing periodic and normal boundary conditions. Boundary normal can be chosen arbitrarily, periodic nodes are only used as value holders.

During this process we also store periodic node indices in $periodicNodes$ and indices of the original nodes in $interiorMapping$.

<syntaxhighlight lang="cpp">
mm::Range<int> periodicNodes;
mm::Range<int> interiorMapping;
mm::Vec2d normal(0, 0); // not used in this case and can be chosen arbitrarily
for (auto idx : domain.interior()) {
const auto& pos = domain.pos(idx);
// col1: distance to - border, col2: distance to + border, rows represent dimensions
auto borderDistance = borders.colwise() - pos;
const auto& isPeriodic = borderDistance.array().abs() < periodicBandWidth * dx;
for (int i0 = 0; i0 < 2; i0++) {
// check for periodicity in 1st dimension (x)
if (isPeriodic(0, i0)) {
auto periodicPos = pos;
periodicPos(0) = borders(0, (i0 + 1) % 2) - borderDistance(0, i0);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
for (int i1 = 0; i1 < 2; i1++) {
// check for periodicity in 2nd dimension (y)
if (i0 == 0 && isPeriodic(1, i1)) {
auto periodicPos = pos;
periodicPos(1) = borders(1, (i1 + 1) % 2) - borderDistance(1, i1);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
// check for periodicity in both dimensions (corners)
if (isPeriodic(0, i0) && isPeriodic(1, i1)) {
auto periodicPos = pos;
periodicPos(0) = borders(0, (i0 + 1) % 2) - borderDistance(0, i0);
periodicPos(1) = borders(1, (i1 + 1) % 2) - borderDistance(1, i1);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
}
}
}
</syntaxhighlight>

== Solution ==
<figure id="fig:timeEvolution">
[[File:TimeEvolution.mp4|thumb|upright=2|<caption>Time evolution of the concentration. Initially mixed fluid separates into domains that grow with time.
</caption>]]
</figure>
We need to handle two uncommon terms to solve the equation. Biharmonic ($\nabla^4$) operator can be implemented with custom operators as shown in [[Customization]] example. Nonlinear ($\nabla^2 c^3$) term proves to be more problematic, as it can not be solved implicitly. We either need to solve the entire equation with explicit time stepping or include the nonlinear term in the right hand side of our system of equations.

Initial concentration is random.

=== Explicit solution ===
Explicit solution is simpler to implement and faster when we want to solve the equation for very small time steps.

We use Euler method for time stepping and copy interior node values to their periodic counterparts at every step.
<syntaxhighlight lang="cpp">
double t = 0;
Eigen::VectorXd bracket;
while (t <= tEnd) {
bracket = concentration.array() * concentration.array() * concentration.array() - concentration.array();
for (int i : allNodes) {
bracket[i] -= L * expOp.lap(concentration, i);
}
for (int i : interiorNodes) {
concentration[i] = concentration[i] + dt * D * expOp.lap(bracket, i);
}
t += dt;

for (int i = 0; i < periodicSize; i++) {
concentration[periodicNodes[i]] = concentration[interiorMapping[i]];
}
}
</syntaxhighlight>

=== Mixed solution ===
Mixed solution allows for much longer (several orders) time steps as we solve the majority of our equation implicitly. Explicit computation of the nonlinear term will always lead to inaccuracies, especially for long $dt$.

We set the equation on interior and identities on boundary.
<syntaxhighlight lang="cpp">
for (int i : interiorNodes) {
impOp.value(i) + D * dt * impOp.lap(i) + D * dt * L * impOp.apply<Biharmonic<dim>>(i) = rhs[i];
}
for (int i = 0; i < periodicSize; i++) {
impOp.value(periodicNodes[i]) = rhs[interiorMapping[i]];
}

Eigen::SparseLU<decltype(M)> solver;
solver.compute(M);
</syntaxhighlight>

We have to explicitly compute nonlinear $\nabla^2 c^3$ term at every step and add it to the rhs. Once we solve the system for new concentration we have to copy interior node values to their periodic counterparts.

<syntaxhighlight lang="cpp">
/// Time stepping.
double t = 0;
Eigen::VectorXd cubedConcentration;
while (t <= tEnd) {
/// Prepare rhs.
cubedConcentration = concentration.array() * concentration.array() * concentration.array();
for (int i : allNodes) {
concentration[i] = D * dt * expOp.lap(cubedConcentration, i) + concentration[i];
}

concentration = solver.solve(concentration);
t += dt;

/// Copy values from interior nodes to their periodic representations.
for (int i = 0; i < periodicSize; i++) {
concentration[periodicNodes[i]] = concentration[interiorMapping[i]];
}
}

</syntaxhighlight>

The complete example is avalible on [https://gitlab.com/e62Lab/medusa/blob/dev/examples/cahnHilliard_equation git].

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

File:TimeEvolution.mp4

2020-02-11T14:10:21Z

MihaR: MihaR uploaded a new version of File:TimeEvolution.mp4

File:TimeEvolution.mp4

2020-02-11T13:42:26Z

MihaR:

Cahn-Hilliard equation

2020-02-11T13:38:02Z

MihaR:

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

In this example we will show how to solve the [https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation Cahn-Hilliard] equation in Medusa. The solution provides an implementation of periodic boundary conditions and shows how to combine implicit and explicit operators. This example uses a custom biharmonic operator, so we recommend firstly reading the tutorial on [[Customization]].

== Cahn-Hilliard equation ==
The Cahn-Hilliard equation describes the separation of binary fluid into pure domains. We define concentration, a scalar field $c(x, y) \in [-1, 1]$ to describe the fluid, with values $c \pm 1$ representing pure phases. In this case the equation can be written as

\[
\frac{\partial c}{\partial t} = D \nabla^2 (c^3 - c - \gamma \nabla^2 c),
\]
where $D$ is the diffusion constant and $\sqrt{\gamma}$ is the characteristic transition length between domains.

We want to avoid any boundary effects and solve the equation for infinite fluid which can be approximated by using periodic boundary conditions.

== Periodic boundary conditions ==

<figure id="fig:periodicBands">
[[File:PeriodicBands.png|thumb|upright=2|<caption>Periodic band display.

Left subplot shows an example of filled domain with $border = 10$ and $dx = 0.2$. Interior nodes and their periodic counterparts share colors to illustrate periodicity.

Right subplot shows the neighbourhood of bottom left corner node $(-border, -border)$, with its support nodes highlighted in red. This example uses same fill parameters as the left figure and a support size of 100.
</caption>]]
</figure>

The simplest method of implementing a periodic boundary is adding a band of periodic nodes that mirror values from nodes that lie close to the opposite boundary. We have to create periodic nodes based on locations of interior nodes and construct a mapping that will allow us to update their values based on their interior counterparts. Width of the periodic band should be chosen with support size in mind. Periodic band should be wide enough to include all nodes required to calculate interior values. An example that can help us determine the required width of the periodic band is shown on the right subplot of <xr id="fig:periodicBands"/> and shows that we need a periodic band width of $6 dx$ for a support size of 100.

We start by creating a square domain centered on 0 and uniformly filling it with nodes. $border$ parameter determines the square size while $dx$, the typical distance between nodes, determines the node density.

<syntaxhighlight lang="cpp">
mm::BoxShape<vec_t> box(-border, border);
mm::DomainDiscretization<vec_t> domain = box.discretizeBoundaryWithStep(dx);
mm::GeneralFill<vec_t> fill;
fill.seed(randomSeed);
domain.fill(fill, dx);
</syntaxhighlight>

Afterwards we remove boundary nodes at positive $border$ and move others into interior. This creates a $[-border, border) \times [-border, border)$ interior that can be made periodic.

<syntaxhighlight lang="cpp">
Eigen::Matrix2d borders; // col1: - border, col2: + border, rows for dimensions
borders << domain.shape().bbox().first, domain.shape().bbox().second;
for (auto idx : domain.boundary()) {
const auto& pos = domain.pos(idx);
if (pos[0] != borders(0, 1) && pos[1] != borders(1, 1)) {
domain.changeToInterior(idx, pos, 1);
}
}
domain.removeBoundaryNodes();
</syntaxhighlight>

Now we are ready to find all nodes that lie sufficiently close to the border and create their periodic counterparts. In general rectangle case we would need to check all dimesions and create periodic nodes in all subsets of dimensions where that node lies close to the border. In this 2D case that simplifies to checking whether node lies close to the border in $x$, $y$ or both, which would mean a periodic corner. We could simplify our code and ignore the small periodic corners but that would lead to slightly worse results.

Periodic boundary nodes are created with type -2 so that we can still differentiate between them and any legitimate boundary nodes we would have in case of mixing periodic and normal boundary conditions. Boundary normal can be chosen arbitrarily, periodic nodes are only used as value holders.

During this process we also store periodic node indices in $periodicNodes$ and indices of the original nodes in $interiorMapping$.

<syntaxhighlight lang="cpp">
mm::Range<int> periodicNodes;
mm::Range<int> interiorMapping;
mm::Vec2d normal(0, 0); // not used in this case and can be chosen arbitrarily
for (auto idx : domain.interior()) {
const auto& pos = domain.pos(idx);
// col1: distance to - border, col2: distance to + border, rows represent dimensions
auto borderDistance = borders.colwise() - pos;
const auto& isPeriodic = borderDistance.array().abs() < periodicBandWidth * dx;
for (int i0 = 0; i0 < 2; i0++) {
// check for periodicity in 1st dimension (x)
if (isPeriodic(0, i0)) {
auto periodicPos = pos;
periodicPos(0) = borders(0, (i0 + 1) % 2) - borderDistance(0, i0);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
for (int i1 = 0; i1 < 2; i1++) {
// check for periodicity in 2nd dimension (y)
if (i0 == 0 && isPeriodic(1, i1)) {
auto periodicPos = pos;
periodicPos(1) = borders(1, (i1 + 1) % 2) - borderDistance(1, i1);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
// check for periodicity in both dimensions (corners)
if (isPeriodic(0, i0) && isPeriodic(1, i1)) {
auto periodicPos = pos;
periodicPos(0) = borders(0, (i0 + 1) % 2) - borderDistance(0, i0);
periodicPos(1) = borders(1, (i1 + 1) % 2) - borderDistance(1, i1);
periodicNodes.push_back(domain.addBoundaryNode(periodicPos, -2, normal));
interiorMapping.push_back(idx);
}
}
}
}
</syntaxhighlight>

== Solution ==
<figure id="fig:timeEvolution">
[[File:TimeEvolution.mp4|thumb|upright=2|<caption>Time evolution of the concentration. Initially mixed fluid separates into domains that grow with time.
</caption>]]
</figure>
We need to handle two uncommon terms to solve the equation. Biharmonic ($\nabla^4$) operator can be implemented with custom operators as shown in [[Customization]] example. Nonlinear ($\nabla^2 c^3$) term proves to be more problematic, as it can not be solved implicitly. We either need to solve the entire equation with explicit time stepping or include the nonlinear term in the right hand side of our system of equations.

Initial concentration is random.

=== Explicit solution ===
Explicit solution is simpler to implement and faster when we want to solve the equation for very small time steps.

We use Euler method for time stepping and copy interior node values to their periodic counterparts at every step.
<syntaxhighlight lang="cpp">
double t = 0;
Eigen::VectorXd bracket;
while (t <= tEnd) {
bracket = concentration.array() * concentration.array() * concentration.array() - concentration.array();
for (int i : allNodes) {
bracket[i] -= L * expOp.lap(concentration, i);
}
for (int i : interiorNodes) {
concentration[i] = concentration[i] + dt * D * expOp.lap(bracket, i);
}
t += dt;

for (int i = 0; i < periodicSize; i++) {
concentration[periodicNodes[i]] = concentration[interiorMapping[i]];
}
}
</syntaxhighlight>

=== Mixed solution ===
Mixed solution allows for much longer (several orders) time steps as we solve the majority of our equation implicitly. Explicit computation of the nonlinear term will always lead to inaccuracies, especially for long $dt$.

We set the equation on interior and identities on boundary.
<syntaxhighlight lang="cpp">
for (int i : interiorNodes) {
impOp.value(i) + D * dt * impOp.lap(i) + D * dt * L * impOp.apply<Biharmonic<dim>>(i) = rhs[i];
}
for (int i = 0; i < periodicSize; i++) {
impOp.value(periodicNodes[i]) = rhs[interiorMapping[i]];
}

Eigen::SparseLU<decltype(M)> solver;
solver.compute(M);
</syntaxhighlight>

We have to explicitly compute nonlinear $\nabla^2 c^3$ term at every step and add it to the rhs. Once we solve the system for new concentration we have to copy interior node values to their periodic counterparts.

<syntaxhighlight lang="cpp">
/// Time stepping.
double t = 0;
Eigen::VectorXd cubedConcentration;
while (t <= tEnd) {
/// Prepare rhs.
cubedConcentration = concentration.array() * concentration.array() * concentration.array();
for (int i : allNodes) {
concentration[i] = D * dt * expOp.lap(cubedConcentration, i) + concentration[i];
}

concentration = solver.solve(concentration);
t += dt;

/// Copy values from interior nodes to their periodic representations.
for (int i = 0; i < periodicSize; i++) {
concentration[periodicNodes[i]] = concentration[interiorMapping[i]];
}
}

</syntaxhighlight>

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

Cahn-Hilliard equation

2020-02-11T00:51:42Z

MihaR:

File:PeriodicBands.png

2020-02-11T00:18:36Z

MihaR:

Cahn-Hilliard equation

2020-02-11T00:17:48Z

MihaR:

Cahn-Hilliard equation

2020-02-07T15:40:23Z

MihaR:

Cahn-Hilliard equation

2020-02-07T13:52:21Z

MihaR:

Cahn-Hilliard equation

2020-02-07T02:09:07Z

MihaR:

Cahn-Hilliard equation

2020-02-07T02:07:55Z

MihaR: MihaR moved page CahnHilliard equation to Cahn-Hilliard equation

CahnHilliard equation

2020-02-07T02:07:55Z

MihaR: MihaR moved page CahnHilliard equation to Cahn-Hilliard equation

#REDIRECT [[Cahn-Hilliard equation]]

Cahn-Hilliard equation

2020-02-07T02:05:51Z

MihaR: Created page with "Go back to Examples. In this example we will show how to solve the [https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation Cahn-Hilliard] equation i..."