http://e6.ijs.si/medusa/wiki/api.php?action=feedcontributions&user=Zigavaupotic&feedformat=atom 2026-04-15T19:39:29Z User contributions MediaWiki 1.27.1 http://e6.ijs.si/medusa/wiki/index.php?title=Customization&diff=3680 2024-07-12T13:31:24Z

<p>Zigavaupotic: /* Custom operators: Hyperviscosity operator */</p> <hr /> <div>Medusa library support users defining custom basis types, weights, operators and more, as long as they conform to the prescribed interfaces, given in the [http://e6.ijs.si/medusa/docs/html/concepts.html Concepts page].<br /> <br /> == Custom stencil selection ==<br /> TODO<br /> <br /> == Custom RBF definition ==<br /> TODO<br /> <br /> == Custom operators: Hyperviscosity operator ==<br /> Medusa library also allows custom linear operators $ \mathcal{L} $ to be used. This functionality is achieved by adding a child struct that inherits the struct Operator. The virtual functions of struct Operator that need to be defined are &lt;code&gt;apply&lt;/code&gt; and &lt;code&gt;applyAt0&lt;/code&gt;, the functions require one to also specify the basis to which the operator is applied. For example if one was to define the hyperviscosity operator denoted as $\nabla^{2\alpha}$, we would say that the RBF-FD approximation to the operator applied to function $u$ is,<br /> $$<br /> \nabla^{2\alpha} u \approx \sum_{i=0}^n u_i \nabla^{2\alpha} \Psi_i (u),<br /> $$<br /> where $\Psi_i (u)$ are cardinal functions and $\alpha$ is the order of hyperviscosity. The evaluated cardinal functions are called weights. The weights are obtained by using a system of RBF coupled with monomial augmentation. Therefore, our operator must have a definition of the function &lt;code&gt;apply&lt;/code&gt; or &lt;code&gt; applyAt0 &lt;/code&gt; for monomial and RBF basis. We define the struct, as<br /> <br /> &lt;syntaxhighlight lang=&quot;cpp&quot;&gt;<br /> template &lt;int dim, int order&gt;<br /> struct Hyperviscosity : public mm::Operator&lt;Hyperviscosity &lt;dim, order&gt;&gt; <br /> &lt;/syntaxhighlight&gt;</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Customization&diff=3679 2024-07-12T13:25:33Z

<p>Zigavaupotic: /* Custom operators: Hyperviscosity operator */</p> <hr /> <div>Medusa library support users defining custom basis types, weights, operators and more, as long as they conform to the prescribed interfaces, given in the [http://e6.ijs.si/medusa/docs/html/concepts.html Concepts page].<br /> <br /> == Custom stencil selection ==<br /> TODO<br /> <br /> == Custom RBF definition ==<br /> TODO<br /> <br /> == Custom operators: Hyperviscosity operator ==<br /> Medusa library also allows custom linear operators $ \mathcal{L} $ to be used. This functionality is achieved by adding a child struct that inherits the struct Operator. The virtual functions of struct Operator that need to be defined are ''apply'' and ''applyAt0'', the functions require one to also specify the basis to which the operator is applied. For example if one was to define the hyperviscosity operator denoted as $\nabla^{2\alpha}$, we would say that the RBF-FD approximation is,<br /> $$<br /> \nabla^{2\alpha} u \approx \sum_{i=0}^n u_i \nabla^{2\alpha} \Psi_i (u),<br /> $$<br /> where $\Psi_i (u)$ are cardinal functions. The evaluated cardinal functions are called weights. The weights are obtained using a system of RBF coupled with monomial augmentation. Therefore, our operator must have a definition of the function ''apply'' for monomial and RBF basis.</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Customization&diff=3678 2024-07-12T12:27:59Z

<p>Zigavaupotic: /* Custom operators: Hyperviscosity operator */</p> <hr /> <div>Medusa library support users defining custom basis types, weights, operators and more, as long as they conform to the prescribed interfaces, given in the [http://e6.ijs.si/medusa/docs/html/concepts.html Concepts page].<br /> <br /> == Custom stencil selection ==<br /> TODO<br /> <br /> == Custom RBF definition ==<br /> TODO<br /> <br /> == Custom operators: Hyperviscosity operator ==<br /> <br /> Medusa library also allows custom linear operators $ \mathcal{L} $ to be used. This functionality is achieved by adding a child struct that inherits the struct Operator. The virtual functions of struct Operator that need to be defined are ''apply'' and ''applyAt0'', the functions require one to also specify the basis to which the operator is applied. For example if one was to define the hyperviscosity operator denoted as $\nabla^{2\alpha}$, we would say that the RBF-FD approximation is,<br /> $$<br /> \nabla^{2\alpha} u \approx \sum_{i=0}^n u_i \nabla^{2\alpha} \Psi_i (u),<br /> $$<br /> where $\Psi_i (u)$ are cardinal functions</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Customization&diff=3677 2024-07-12T12:27:06Z

<p>Zigavaupotic: /* Custom operators: Hyperviscosity operator */</p> <hr /> <div>Medusa library support users defining custom basis types, weights, operators and more, as long as they conform to the prescribed interfaces, given in the [http://e6.ijs.si/medusa/docs/html/concepts.html Concepts page].<br /> <br /> == Custom stencil selection ==<br /> TODO<br /> <br /> == Custom RBF definition ==<br /> TODO<br /> <br /> == Custom operators: Hyperviscosity operator ==<br /> <br /> Medusa library also allows custom linear operators $ \mathcal{L} $ to be used. This functionality is achieved by adding a child struct that inherits the struct Operator. The virtual functions of struct Operator that need to be defined are ''apply'' and ''applyAt0'', the functions require one to also specify the basis to which the operator is applied. For example if one was to define the hyperviscosity operator denoted as $\nabla^{2\alpha}$, we would say that the RBF-FD approximation is,<br /> $$<br /> \nabla^{2\alpha} u \approx \sum_{i=0}^n u_i \nabla^{2\alpha} \Psi_i (u) + \sum_{j=0}^\ell \nabla^{2\alpha}\beta_{i} u_i,<br /> $$<br /> where $\Psi_i (u)$ are cardinal functions and $\beta_i(u_i)$ are the monomial weights</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Customization&diff=3676 2024-07-12T12:26:56Z

<p>Zigavaupotic: /* Custom operators: Hyperviscosity operator */</p> <hr /> <div>Medusa library support users defining custom basis types, weights, operators and more, as long as they conform to the prescribed interfaces, given in the [http://e6.ijs.si/medusa/docs/html/concepts.html Concepts page].<br /> <br /> == Custom stencil selection ==<br /> TODO<br /> <br /> == Custom RBF definition ==<br /> TODO<br /> <br /> == Custom operators: Hyperviscosity operator ==<br /> <br /> Medusa library also allows custom linear operators $ \mathcal{L} $ to be used. This functionality is achieved by adding a child struct that inherits the struct Operator. The virtual functions of struct Operator that need to be defined are ''apply'' and ''applyAt0'', the functions require one to also specify the basis to which the operator is applied. For example if one was to define the hyperviscosity operator denoted as $\nabla^{2\alpha}$, we would say that the RBF-FD approximation is,<br /> $$<br /> \nabla^{2\alpha} u \approx \sum_{i=0}^n u_i \nabla^{2\alpha} \Psi_i (u) + \sum_{j=0}^\ell \nabla^{2\alpha}\beta_{i} u_i,<br /> $$<br /> where $\Psi_i (u)$ are cardinal functions and $\beta$ are the monomial weights</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Customization&diff=3675 2024-07-12T12:26:23Z

<p>Zigavaupotic: /* Custom operators: Biharmonic equation */</p> <hr /> <div>Medusa library support users defining custom basis types, weights, operators and more, as long as they conform to the prescribed interfaces, given in the [http://e6.ijs.si/medusa/docs/html/concepts.html Concepts page].<br /> <br /> == Custom stencil selection ==<br /> TODO<br /> <br /> == Custom RBF definition ==<br /> TODO<br /> <br /> == Custom operators: Hyperviscosity operator ==<br /> <br /> Medusa library also allows custom linear operators $ \mathcal{L} $ to be used. This functionality is achieved by adding a child struct that inherits the struct Operator. The virtual functions of struct Operator that need to be defined are ''apply'' and ''applyAt0'', the functions require one to also specify the basis to which the operator is applied. For example if one was to define the hyperviscosity operator denoted as $\nabla^{2\alpha}$, we would say that the RBF-FD approximation is,<br /> $$<br /> \nabla^{2\alpha} u \approx \sum_{i=0}^n u_i \nabla^{2\alpha} \Psi_i (u) + \sum_{j=0}^\ell \nabla^{2\alpha}\beta_{i} u_i<br /> $$</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Hyperviscosity&diff=3672 2024-03-23T14:24:02Z

<p>Zigavaupotic: </p> <hr /> <div>One of the most common drawbacks of the RBF-FD method is its stability. Furthermore, the RBF-FD method, due to scattered nodes, sometimes produces system (differential) matrices with spurious eigenvalues. Those are far more common than with for example FVM or FEM method. <br /> <br /> For such unstable systems, one usually considers numerical diffusion in the form of hyperviscosity. The hyperviscosity scheme introduces higher-order Laplacian operator (i.e. hyperviscosity operator) to the right hand side of the PDE<br /> <br /> $$<br /> +\gamma \Delta^\alpha u,<br /> $$<br /> <br /> where $\gamma$ is a hyperviscosity constant and $\alpha$ is the order of hyperviscosity.</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Hyperviscosity&diff=3671 2024-03-23T14:19:29Z

<p>Zigavaupotic: Created page with &quot;One of the most commonly mentioned drawbacks of the RBF-FD method is its stability. Furthermore, the RBF-FD method, due to scattered nodes, sometimes produces system (differen...&quot;</p> <hr /> <div>One of the most commonly mentioned drawbacks of the RBF-FD method is its stability. Furthermore, the RBF-FD method, due to scattered nodes, sometimes produces system (differential) matrices with spurious eigenvalues. Those are far more common than with for example FVM or FEM method. The instabilities arise fastly even with linear PDEs.<br /> <br /> With meshless methods (especially RBF-FD) one usually stabilises such schemes by adding a higher-order Laplacian (hyperviscosity) operator on the right-hand side.<br /> <br /> $$<br /> +\gamma \Delta^\alpha u<br /> $$<br /> <br /> The final scheme is referred to as hyperviscosity scheme</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Medusa&diff=3569 2024-01-13T13:18:56Z

<p>Zigavaupotic: </p> <hr /> <div>&lt;!--__NOTITLE__--&gt;<br /> '''Welcome to the Medusa wiki. To visit the main website, go to [http://e6.ijs.si/medusa/ http://e6.ijs.si/medusa/].'''<br /> <br /> 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.<br /> We call this open source meshless project [http://e6.ijs.si/medusa/ Medusa: Coordinate Free Meshless Method implementation (MM)].<br /> <br /> 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]]<br /> <br /> This wiki site is meant for more relaxed discussions about general principles, possible and already implemented applications, preliminary analyses, etc.<br /> 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 :). <br /> <br /> == Documentation ==<br /> * [https://gitlab.com/e62Lab/medusa Code on Gitlab]<br /> * [[How to build | Installation and building]]<br /> * [[Including this library in your project | Including this library in your project]]<br /> * [[Testing | Running tests]]<br /> * [http://e6.ijs.si/medusa/docs/ Technical documentation]<br /> * [[Coding style | Coding style]]<br /> * [[Wiki editing guide | Wiki editing and backup guide]]<br /> <br /> == Building blocks ==<br /> 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. <br /> * [[Positioning of computational nodes]] <br /> * [[Relaxation of the nodal distribution]]<br /> * [[Refinement of the nodal distribution]]<br /> * [[k-d tree|''k''-d tree and other spatial search structures]] <br /> * [[Solving system | Solving linear system - including over and underdetermined systems]]<br /> * [[Weighted Least Squares (WLS)]]<br /> * [[Computation of shape functions]]<br /> * [[Radial basis function-generated finite differences (RBF-FD)]]<br /> * [[Ghost nodes (theory)]]<br /> * [[Integrators for time stepping]]<br /> <br /> == Examples ==<br /> 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.<br /> * [[Philosophy of examples and how to run them]]<br /> * [[Poisson's equation]]<br /> * [[Heat equation]]<br /> * [[Linear elasticity]]<br /> * [[Complex-valued problems]]<br /> * [[Coupled domains]]<br /> * [[Parametric domains | Parametric domains &amp;ndash; Curved surface with variable density]]<br /> * [[NURBS domains | Domains modeled with non-uniform rational basis splines (NURBS)]]<br /> * [[Determining the interior of the domain by oversampling the boundary]]<br /> * [[Computer-aided design - Importing IGES and STEP files]]<br /> * [[Realistic 3D models|Working with 3D surface mesh models ]]<br /> * [[customization | Operator customization]]<br /> * [[Ghost nodes]]<br /> * [[Electromagnetic scattering]]<br /> * [[Schrödinger equation]]<br /> * [[Wave equation]]<br /> * [[Cahn-Hilliard equation]]<br /> * [[Fluid mechanics]]<br /> * [[Solid Mechanics | Solid mechanics]]<br /> <br /> == Discussions / Applications ==<br /> 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. <br /> * [[Convection Diffusion equation | Convection Diffusion equation]]<br /> * [[Adaptivity|H-adaptivity]]<br /> * [[Hp-adaptivity]]<br /> * [[Solid Mechanics]]<br /> ** [[Point contact]]<br /> ** [[Hertzian contact]]<br /> ** [[Cantilever beam]]<br /> ** [[Fretting fatigue case]]<br /> ** [[Plasticity]]<br /> * [[Fluid Mechanics]]<br /> ** [[Lid driven cavity]]<br /> ** [[Burger's equation]]<br /> ** [[de Vahl Davis natural convection test]]<br /> ** [[Natural convection in 3D irregular domain]]<br /> ** [[Natural convection from heated cylinder]]<br /> ** [[Natural convection between concentric cylinders]]<br /> ** [[Non-Newtonian fluid]]<br /> * [[Computational electromagnetics]]<br /> ** [[Triple dielectric step in 1D]]<br /> ** [[Scattering from an infinite cylinder]]<br /> ** [[Point source near an anisotropic lens]]<br /> * Other applications<br /> ** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]<br /> ** [[Heart rate variability detection]]<br /> ** [[Bioheat equation]]<br /> <br /> == Performance analyses ==<br /> * [[Execution on Intel® Xeon Phi™ co-processor]]<br /> * [[1D MLSM and FDM comparison]]<br /> * [[:File:tech_report.pdf|Execution overheads due to clumsy types::technical report]] [[File:pdf-file.gif]]<br /> * [[Solving sparse systems]]<br /> * [[Eigen paralelization]]<br /> <br /> == Last changes ==<br /> &lt;news unique=1 limit = 5&gt;<br /> *{{{timeanddate}}} :: {{{title}}} <br /> <br /> &lt;/news&gt;<br /> <br /> == Miscellaneous ==<br /> * FAQ - [[Frequently asked questions]]. <br /> * [[List of wiki contributors]]<br /> * List of library contributors: [http://e6.ijs.si/medusa/about#about-contributors See the official website]<br /> <br /> == References ==<br /> <br /> For all related papers including conference contributions, monographs and book chapters check https://e6.ijs.si/ParallelAndDistributedSystems/publications/<br /> <br /> {{Box-round|title= Selected papers |<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/98533123.pdf M. Depolli, J. Slak, G. Kosec; Parallel domain discretization algorithm for RBF-FD and other meshless numerical methods for solving PDEs, Computers &amp; Structures, 2022 [DOI: 10.1016/j.compstruc.2022.106773]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/69777155.pdf J. Slak, G. Kosec; Medusa : A C++ library for solving PDEs using strong form mesh-free methods, ACM transactions on mathematical software, vol. 47, 2021 [DOI: 10.1145/3450966]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/56730115.pdf U. Duh, G. Kosec, J. Slak; Fast variable density node generation on parametric surfaces with application to mesh-free methods, SIAM journal on scientific computing, vol. 43, 2021 [DOI: 10.1137/20M1325642]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/52715011.pdf M. Jančič, J. Slak, G. Kosec; Monomial augmentation guidelines for RBF-FD from accuracy versus computational time perspective, Journal of scientific computing, vol. 87, 2021 [DOI: 10.1007/s10915-020-01401-y]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32782887.pdf J. Slak, G. Kosec; On generation of node distributions for meshless PDE discretizations, SIAM journal on scientific computing, vol. 41, 2019 [DOI: 10.1137/18M1231456]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32424999.pdf G. Kosec, J. Slak, M. Depolli, R. Trobec, K. Pereira, S. Tomar, T. Jacquemin, S. Bordas, M. Wahab; Weak and strong from meshless methods for linear elastic problem under fretting contact conditions, Tribology international, vol. 138, 2019]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32230439.pdf J. Slak, G. Kosec; Adaptive radial basis function-generated finite differences method for contact problems, International journal for numerical methods in engineering, vol. 119, 2019 [DOI: 10.1002/nme.6067]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32388135.pdf M. Maksić, V. Djurica, A. Souvent, J. Slak, M. Depolli, G. Kosec; Cooling of overhead power lines due to the natural convection, International journal of electrical power &amp; energy systems, 2019]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/31107623.pdf J. Slak, G. Kosec; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain, Engineering analysis with boundary elements, vol. 100, 2019]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/29639719.pdf M. Depolli, G. Kosec; 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]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/29512743.pdf G. Kosec; A local numerical solution of a fluid-flow problem on an irregular domain, Advances in engineering software, vol. 120, 2018 [DOI: 10.1016/j.advengsoft.2016.05.010]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/27912487.pdf G. Kosec, R. Trobec; Simulation of semiconductor devices with a local numerical approach, Engineering analysis with boundary elements, 2015 [DOI: 10.1016/j.enganabound.2014.07.013]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/3218939.pdf G. Kosec, B. Šarler; Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method, Engineering analysis with boundary elements]<br /> <br /> [http://comms.ijs.si/~gkosec/data/papers/27339815.pdf G. Kosec, M. Depolli, A. Rashkovska, R. Trobec; Super linear speedup in a local parallel meshless solution of thermo-fluid problem, Computers &amp; Structures, vol. 133, 2014]<br /> <br /> [http://comms.ijs.si/~gkosec/data/papers/26785063.pdf G. Kosec, P. Zinterhof; Local strong form meshless method on multiple Graphics Processing Units, Computer modeling in engineering &amp; sciences, vol. 91, 2013]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/2599419.pdf G. Kosec, B. Šarler; Solution of a low Prandtl number natural convection benchmark by a local meshless method, International journal of numerical methods for heat &amp; fluid flow]<br /> <br /> [http://comms.ijs.si/~gkosec/data/papers/EABE2499.pdf R. Trobec, G. Kosec, M. Šterk, B. Šarler; Comparison of local weak and strong form meshless methods for 2-D diffusion equation, Engineering analysis with boundary elements, vol. 36, 2012]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/1905659.pdf G. Kosec, M. Založnik, B. Šarler, H. Combeau; A meshless approach towards solution of macrosegregation phenomena, Computers, materials &amp; continua : CMC, vol. 22, 2011 ]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/992507.pdf G. Kosec, B. Šarler; Solution of thermo-fluid problems by collocation with local pressure correction, International journal of numerical methods for heat &amp; fluid flow, vol.18, 2008]<br /> <br /> R. Trobec, G. Kosec; Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods, 2015<br /> <br /> }}<br /> <br /> ==Related pages==<br /> * http://e6.ijs.si/ParallelAndDistributedSystems/products/medusa/</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Medusa&diff=3568 2024-01-13T13:18:35Z

<p>Zigavaupotic: </p> <hr /> <div>&lt;!--__NOTITLE__--&gt;<br /> '''Welcome to the Medusa wiki. To visit the main website, go to [http://e6.ijs.si/medusa/ http://e6.ijs.si/medusa/].'''<br /> <br /> 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.<br /> We call this open source meshless project [http://e6.ijs.si/medusa/ Medusa: Coordinate Free Meshless Method implementation (MM)].<br /> <br /> 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]]<br /> <br /> This wiki site is meant for more relaxed discussions about general principles, possible and already implemented applications, preliminary analyses, etc.<br /> 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 :). <br /> <br /> == Documentation ==<br /> * [https://gitlab.com/e62Lab/medusa Code on Gitlab]<br /> * [[How to build | Installation and building]]<br /> * [[Including this library in your project | Including this library in your project]]<br /> * [[Testing | Running tests]]<br /> * [http://e6.ijs.si/medusa/docs/ Technical documentation]<br /> * [[Coding style | Coding style]]<br /> * [[Wiki editing guide | Wiki editing and backup guide]]<br /> <br /> == Building blocks ==<br /> 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. <br /> * [[Positioning of computational nodes]] <br /> * [[Relaxation of the nodal distribution]]<br /> * [[Refinement of the nodal distribution]]<br /> * [[k-d tree|''k''-d tree and other spatial search structures]] <br /> * [[Solving system | Solving linear system - including over and underdetermined systems]]<br /> * [[Weighted Least Squares (WLS)]]<br /> * [[Computation of shape functions]]<br /> * [[Radial basis function-generated finite differences (RBF-FD)]]<br /> * [[Ghost nodes (theory)]]<br /> * [[Integrators for time stepping]]<br /> <br /> == Examples ==<br /> 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.<br /> * [[Philosophy of examples and how to run them]]<br /> * [[Poisson's equation]]<br /> * [[Heat equation]]<br /> * [[Linear elasticity]]<br /> * [[Complex-valued problems]]<br /> * [[Coupled domains]]<br /> * [[Parametric domains | Parametric domains &amp;ndash; Curved surface with variable density]]<br /> * [[NURBS domains | Domains modeled with non-uniform rational basis splines (NURBS)]]<br /> * [[Determining the interior of the domain by oversampling the boundary]]<br /> * [[Computer-aided design - Importing IGES and STEP files]]<br /> * [[Realistic 3D models|Working with 3D surface mesh models ]]<br /> * [[customization | Operator customization]]<br /> * [[Ghost nodes]]<br /> * [[Electromagnetic scattering]]<br /> * [[Schrödinger equation]]<br /> * [[Wave equation]]<br /> * [[Cahn-Hilliard equation]]<br /> * [[Fluid mechanics]]<br /> * [[Solid Mechanics | Solid mechanics]]<br /> <br /> == Discussions / Applications ==<br /> 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. <br /> * [[Convection Diffusion equation | Convection Diffusion equation]]<br /> * [[Adaptivity|H-adaptivity]]<br /> * [[Hp-adaptivity]]<br /> * [[Solid Mechanics]]<br /> ** [[Point contact]]<br /> ** [[Hertzian contact]]<br /> ** [[Cantilever beam]]<br /> ** [[Fretting fatigue case]]<br /> ** [[Plasticity]]<br /> * [[Fluid Mechanics]]<br /> ** [[Lid driven cavity]]<br /> ** [[Burger's Equation]]<br /> ** [[de Vahl Davis natural convection test]]<br /> ** [[Natural convection in 3D irregular domain]]<br /> ** [[Natural convection from heated cylinder]]<br /> ** [[Natural convection between concentric cylinders]]<br /> ** [[Non-Newtonian fluid]]<br /> * [[Computational electromagnetics]]<br /> ** [[Triple dielectric step in 1D]]<br /> ** [[Scattering from an infinite cylinder]]<br /> ** [[Point source near an anisotropic lens]]<br /> * Other applications<br /> ** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]<br /> ** [[Heart rate variability detection]]<br /> ** [[Bioheat equation]]<br /> <br /> == Performance analyses ==<br /> * [[Execution on Intel® Xeon Phi™ co-processor]]<br /> * [[1D MLSM and FDM comparison]]<br /> * [[:File:tech_report.pdf|Execution overheads due to clumsy types::technical report]] [[File:pdf-file.gif]]<br /> * [[Solving sparse systems]]<br /> * [[Eigen paralelization]]<br /> <br /> == Last changes ==<br /> &lt;news unique=1 limit = 5&gt;<br /> *{{{timeanddate}}} :: {{{title}}} <br /> <br /> &lt;/news&gt;<br /> <br /> == Miscellaneous ==<br /> * FAQ - [[Frequently asked questions]]. <br /> * [[List of wiki contributors]]<br /> * List of library contributors: [http://e6.ijs.si/medusa/about#about-contributors See the official website]<br /> <br /> == References ==<br /> <br /> For all related papers including conference contributions, monographs and book chapters check https://e6.ijs.si/ParallelAndDistributedSystems/publications/<br /> <br /> {{Box-round|title= Selected papers |<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/98533123.pdf M. Depolli, J. Slak, G. Kosec; Parallel domain discretization algorithm for RBF-FD and other meshless numerical methods for solving PDEs, Computers &amp; Structures, 2022 [DOI: 10.1016/j.compstruc.2022.106773]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/69777155.pdf J. Slak, G. Kosec; Medusa : A C++ library for solving PDEs using strong form mesh-free methods, ACM transactions on mathematical software, vol. 47, 2021 [DOI: 10.1145/3450966]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/56730115.pdf U. Duh, G. Kosec, J. Slak; Fast variable density node generation on parametric surfaces with application to mesh-free methods, SIAM journal on scientific computing, vol. 43, 2021 [DOI: 10.1137/20M1325642]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/52715011.pdf M. Jančič, J. Slak, G. Kosec; Monomial augmentation guidelines for RBF-FD from accuracy versus computational time perspective, Journal of scientific computing, vol. 87, 2021 [DOI: 10.1007/s10915-020-01401-y]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32782887.pdf J. Slak, G. Kosec; On generation of node distributions for meshless PDE discretizations, SIAM journal on scientific computing, vol. 41, 2019 [DOI: 10.1137/18M1231456]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32424999.pdf G. Kosec, J. Slak, M. Depolli, R. Trobec, K. Pereira, S. Tomar, T. Jacquemin, S. Bordas, M. Wahab; Weak and strong from meshless methods for linear elastic problem under fretting contact conditions, Tribology international, vol. 138, 2019]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32230439.pdf J. Slak, G. Kosec; Adaptive radial basis function-generated finite differences method for contact problems, International journal for numerical methods in engineering, vol. 119, 2019 [DOI: 10.1002/nme.6067]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/32388135.pdf M. Maksić, V. Djurica, A. Souvent, J. Slak, M. Depolli, G. Kosec; Cooling of overhead power lines due to the natural convection, International journal of electrical power &amp; energy systems, 2019]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/31107623.pdf J. Slak, G. Kosec; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain, Engineering analysis with boundary elements, vol. 100, 2019]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/29639719.pdf M. Depolli, G. Kosec; 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]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/29512743.pdf G. Kosec; A local numerical solution of a fluid-flow problem on an irregular domain, Advances in engineering software, vol. 120, 2018 [DOI: 10.1016/j.advengsoft.2016.05.010]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/27912487.pdf G. Kosec, R. Trobec; Simulation of semiconductor devices with a local numerical approach, Engineering analysis with boundary elements, 2015 [DOI: 10.1016/j.enganabound.2014.07.013]]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/3218939.pdf G. Kosec, B. Šarler; Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method, Engineering analysis with boundary elements]<br /> <br /> [http://comms.ijs.si/~gkosec/data/papers/27339815.pdf G. Kosec, M. Depolli, A. Rashkovska, R. Trobec; Super linear speedup in a local parallel meshless solution of thermo-fluid problem, Computers &amp; Structures, vol. 133, 2014]<br /> <br /> [http://comms.ijs.si/~gkosec/data/papers/26785063.pdf G. Kosec, P. Zinterhof; Local strong form meshless method on multiple Graphics Processing Units, Computer modeling in engineering &amp; sciences, vol. 91, 2013]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/2599419.pdf G. Kosec, B. Šarler; Solution of a low Prandtl number natural convection benchmark by a local meshless method, International journal of numerical methods for heat &amp; fluid flow]<br /> <br /> [http://comms.ijs.si/~gkosec/data/papers/EABE2499.pdf R. Trobec, G. Kosec, M. Šterk, B. Šarler; Comparison of local weak and strong form meshless methods for 2-D diffusion equation, Engineering analysis with boundary elements, vol. 36, 2012]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/1905659.pdf G. Kosec, M. Založnik, B. Šarler, H. Combeau; A meshless approach towards solution of macrosegregation phenomena, Computers, materials &amp; continua : CMC, vol. 22, 2011 ]<br /> <br /> [https://e6.ijs.si/ParallelAndDistributedSystems/publications/992507.pdf G. Kosec, B. Šarler; Solution of thermo-fluid problems by collocation with local pressure correction, International journal of numerical methods for heat &amp; fluid flow, vol.18, 2008]<br /> <br /> R. Trobec, G. Kosec; Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods, 2015<br /> <br /> }}<br /> <br /> ==Related pages==<br /> * http://e6.ijs.si/ParallelAndDistributedSystems/products/medusa/</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Burgers%27_equation&diff=3567 2024-01-13T12:32:10Z

<p>Zigavaupotic: </p> <hr /> <div>Let us consider the Burger's equation, which describes the dynamics of viscous fluid without the effects of pressure. Despite being unrealistic it is the simplest description of advective flow with diffusive effects of viscosity. It has the following form<br /> \begin{equation}<br /> \frac{\partial \b{u}}{\partial t} + (\b{u}\cdot\nabla)\b{u} = \nu \nabla^2 \b{u}<br /> \end{equation}<br /> <br /> where $\b{u}$ is a velocity field. For time stepping the advection term must be linearized. There are several tactics of linearization, but we choose the simplest one, which is using the current velocity as an approximant of the velociy in the next time step. Using the implicit Euler method we arrive at the equation<br /> <br /> \begin{equation}<br /> \b{u}^{n+1}_i + [\b{u}^{n}_i\cdot\nabla_i\b{u}^{n+1} - \nu\nabla^2_i \b{u}^{n+1}]\Delta t = \b{u}^{n}_i<br /> \end{equation}<br /> <br /> where the subscript $i$ is the index of the domain node, superscripts $n$ and $n+1$ denote the current and the next time step respectively, $\Delta t$ is the length of the time step. The notation $\nabla_i$ means the gradient operator at node $i$, not the component of the gradient. Same holds for $\nabla^2_i$. Using the Medusa library, this equation can be solved for $\b{u}^{n+1}_i$.<br /> <br /> =Solution in 1D=<br /> Let's look at an example in 1D with Dirichlet boundary conditions. We are solving the equivalent of equation (2) but with only one spatial variable<br /> $$u^{n+1}_i + [u^{n}_i\frac{\partial u^{n+1}}{\partial x} - \nu\frac{\partial^2 u^{n+1}}{\partial x^2}]\Delta t = u^{n}_i$$<br /> $$u_i = 0 \text{ on } \partial \Omega$$<br /> <br /> This will be rewritten into a linear system of equations with Medusa and then solved with Eigen. Let's look at the code(find the full version in our repository)<br /> <br /> &lt;syntaxhighlight lang=&quot;c++&quot; line&gt;<br /> // Build 1D domain<br /> BoxShape&lt;Vec1d&gt; domain((-l), (l));<br /> auto discretization = domain.discretizeWithStep(dx);<br /> int domain_size = discretization.size();<br /> // Find support nodes<br /> FindClosest find_support(5);<br /> discretization.findSupport(find_support);<br /> &lt;/syntaxhighlight&gt;<br /> <br /> First we construct our 1D domain with lenght 2l, discretize it with a given step dx and find the closest 5 nodes of each node which are called stencil nodes.<br /> <br /> &lt;syntaxhighlight lang=&quot;c++&quot; line&gt;<br /> Polyharmonic&lt;double, 3&gt; p;<br /> Monomials&lt;Vec1d&gt; mon(4);<br /> RBFFD&lt;Polyharmonic&lt;double, 3&gt;, Vec1d, ScaleToClosest&gt;<br /> approx(p, mon);<br /> auto storage = discretization.computeShapes(approx); <br /> VectorXd E1 = VectorXd::Zero(domain_size); // vector for containing the current state<br /> &lt;/syntaxhighlight&gt;<br /> <br /> We will use Polyharmonic Radial Basis Functions(RBFs) of order 3 augmented with monomials up to 3th order for our approximation engine. We use it in conjunction with stencil nodes to compute the stencil weights(also called shape functions) which are used in approximations of differential operators on our domain. Next we set the initial condition to $$u(x, 0) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)}$$ where $a, b, k$ are constants. We chose this function because it has an analytical solution which we will compare our solution to. It is derived using the Cole-Hopf transformation of a solution of the diffusion equation with the initial profile $f(x) = b + a\cos(kx)$ for $b&gt;a$.<br /> <br /> &lt;syntaxhighlight lang=&quot;c++&quot; line&gt;<br /> //Set initial condition <br /> for (int i : discretization.all()){ <br /> E1(i) = 2*visc*a*k*std::exp(-visc*pow(k,2)*t_0)*std::sin(k*discretization.pos(i,0))/(b+a*std::exp(- <br /> visc*pow(k,2)*t_0)*std::cos(k*discretization.pos(i,0)));<br /> &lt;/syntaxhighlight&gt;<br /> <br /> Inside our time-stepping loop, we instantiate the matrix $M$ and the vector $rhs$ representing the linear system. We then construct the differential operators for implicit solving.<br /> &lt;syntaxhighlight lang=&quot;c++&quot; line&gt;<br /> for (tt = 1; tt &lt;= t_steps; ++tt) {<br /> SparseMatrix&lt;double&gt; M(domain_size, domain_size); // construct matrix of apropreate size<br /> VectorXd rhs = VectorXd::Zero(domain_size); // set empty vector for rhs<br /> M.reserve(Range&lt;int&gt;(domain_size, n));<br /> auto op = storage.implicitOperators(M, rhs); // compute the implicit operators <br /> &lt;/syntaxhighlight&gt;<br /> <br /> Now comes the crucial part - setting the equation. In essence, we rewrite the equation (2) in Medusa's terms. Terms with the prefix op. are matrices that represent the appropriate differential operators(op.value is simply the identity matrix).</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Burgers%27_equation&diff=3564 2023-12-27T12:34:28Z

<p>Zigavaupotic: </p> <hr /> <div>Burgers equation<br /> <br /> <br /> <br /> $$\frac{\partial u}{\partial t} + u(\nabla \cdot u) = \nu \nabla^2 u$$</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Burgers%27_equation&diff=3563 2023-12-27T12:34:17Z

<p>Zigavaupotic: </p> <hr /> <div>Burgers equation<br /> <br /> <br /> <br /> $$\frac{\partial u}{\partial t} + u(\nabla \cdot u) = \ni \nabla^2 u$$</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Fluid_Mechanics&diff=3562 2023-12-27T12:32:07Z

<p>Zigavaupotic: </p> <hr /> <div>= Introduction =<br /> Computational fluid dynamics (CFD) is a field of a great interest among researchers in many fields of science, e.g. studying mathematical fundaments of numerical methods, developing novel physical models, improving computer implementations, and many others. Pushing the limits of all involved fields of science helps community to deepen the understanding of several natural and technological phenomena. Weather forecast, ocean dynamics, water transport, casting, various energetic studies, etc., are just few examples where fluid dynamics plays a crucial role. The core problem of the CFD is solving the Navier-Stokes Equation or its variants, e.g. Darcy or Brinkman equation for flow in porous media. Here, we discuss basic algorithms for solving CFD problems. Check reference list on the [[Medusa|main page]] for more details about related work.<br /> <br /> Long story short, we want to solve<br /> \begin{equation}<br /> \frac{\partial \b{v}}{\partial t}+\nabla \cdot ( \rho \b{vv}) = -\nabla p+ \nabla(\mu \nabla\b{v})+\b{f}<br /> \label{NavierStokes}<br /> \end{equation}<br /> also known as a Navier-Stokes equation. <br /> <br /> Note that the $\b{v}\b{v}$ stands for the tensor or dyadic product \[ \b{v}\b{v} = \b{v}\otimes\b{v} = \b{v}\b{v}^\T = \left[ \begin{matrix}<br /> {{v}_{1}}{{v}_{1}} &amp; \cdots &amp; {{v}_{1}}{{v}_{n}} \\<br /> \vdots &amp; \ddots &amp; \vdots \\<br /> {{v}_{n}}{{v}_{1}} &amp; \cdots &amp; {{v}_{n}}{{v}_{n}} \\<br /> \end{matrix} \right]\]<br /> <br /> Combining equation \ref{NavierStokes} and mass contiunity equation<br /> \begin{equation}<br /> \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\bf{u}}) = 0.<br /> \end{equation}<br /> <br /> simplifies advection term into (see [https://en.wikipedia.org/wiki/Derivation_of_the_Navier%E2%80%93Stokes_equations] for derivation)<br /> <br /> \[\frac{\partial \left( \rho \b{v} \right)}{\partial t}+\nabla \cdot \left( \rho \b{vv} \right)=\frac{\partial \left( \rho \b{v} \right)}{\partial t}+(\rho \b{v}\cdot \nabla )\cdot \b{v}. \]<br /> <br /> An example of advection term in 2D would therefore be<br /> \[\left( \b{v}\cdot \nabla \right)\b{v}=\left( \left( \begin{matrix}<br /> u \\<br /> v \\<br /> \end{matrix} \right) \cdot \left( \begin{matrix}<br /> \frac{\partial }{\partial x} \\<br /> \frac{\partial }{\partial y} \\<br /> \end{matrix} \right) \right)\left( \begin{matrix}<br /> u \\<br /> v \\<br /> \end{matrix} \right)=\left( \begin{matrix}<br /> u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} \\<br /> u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} \\<br /> \end{matrix} \right)\]<br /> <br /> <br /> In many cases we are interested in the incompressible fluids (Ma&lt;0.3), reducing the continuity equation to<br /> \begin{equation}<br /> \nabla \cdot \b{v}=0<br /> \label{contuinity}<br /> \end{equation}<br /> <br /> <br /> The goal of CFD is to solve system \ref{NavierStokes} and \ref{contuinity}. It is obvious that a special treatment will be needed to couple both equations. In the following discussion we cover some basic approaches, how this can be accomplished.<br /> <br /> = Solutions algorithms =<br /> == Artificial compressibility method ==<br /> The simplest, completely explicit approach, is an artificial compressibility method (ACM), where a compressibility term is included in the mass continuity<br /> \[\frac{\partial \b{v}}{\partial t}+(\b{v}\cdot\nabla )\b{v}=-\frac{1}{\rho }\nabla p+\nu {{\nabla }^{2}}\b{v}+\b{f}\]<br /> \[\frac{ 1 }{ \rho } \frac{\partial \rho }{\partial t}+\nabla \cdot \b{v}=0\]<br /> \[\frac{ 1 }{ \rho } \frac{\partial \rho }{\partial p}\frac{\partial p}{\partial t}+\nabla \cdot \b{v}=0\]<br /> Now, the above system can be solved directly.<br /> <br /> The addition of the time derivative of the pressure term physically means that waves of finite speed (the propagation of which depends on the magnitude of the ACM)<br /> are introduced into the flow field as a mean to distribute the pressure within the domain. In a true<br /> incompressible flow, the pressure field is affected instantaneously throughout the whole domain. In ACM there is a time delay between the flow disturbance and its effect on the<br /> pressure field. Upon rearranging the equation yields<br /> \[\frac{\partial p}{\partial t}+\rho {{C}^{2}}\nabla \cdot \b{v}=0\]<br /> where the continuity equation is perturbed by the quantity $\frac{\partial p}{\partial t}$ denominated herein<br /> as the AC parameter/artificial sound speed recognized by<br /> $C$ [m/s] - speed of sound<br /> \[\frac{1}{C^2}=\frac{\partial \rho }{\partial p}\]<br /> Or in other words<br /> \[C^2=\left( \frac{\partial p}{\partial \rho}\right)_S\]<br /> where $\rho$ is the density of the material. It follows, by replacing partial derivatives, that the isentropic compressibility can be expressed as:<br /> \[\beta =\frac{1}{\rho {{C}^{2}}}\]<br /> The evaluation of the local ACM parameter in incompressible flows is inspired by the<br /> speed of sound computations in compressible flows (for instance, from the perfect gas law).<br /> However, in the incompressible flow situation, employing such a relation is difficult, but an artificial<br /> relation can be developed from the convective and diffusive velocities.<br /> Reverting to the justification of continuity modification, it can be immediately seen that the<br /> artificial sound speed must be sufficiently large to have a significant regularizing effect and at<br /> the same time must be as small as possible to minimizing perturbations on the incompressibility<br /> equation. Therefore, $C$ influences the convergence rate and stability of the solution method. In other words,<br /> assists in reducing large disparity in the eigenvalues, leading to a well-conditioned system. <br /> The $C$ can be '''estimated''' with<br /> <br /> \[ C = \beta \max \left( \left|\b{v}\right|_2, \left|\b{v}_{ref}\right|_2 \right),\]<br /> where $\b{v}_{ref}$ stands for a reference velocity. <br /> Values for $\beta$ in the range of 1–10 are recommended for better convergence to the steady state at which the<br /> mass conservation is enforced. In addition, Equation ensures that $C$ does not reach zero at stagnation points<br /> that cause instabilities in pseudo-time, effecting convergence<br /> <br /> Note, that for more complex simulation an internal iteration loops is required before marching in time. <br /> <br /> <br /> &lt;figure id=&quot;ACM&quot;&gt;<br /> [[File:ACM BlockDiagram.png|600px]]<br /> &lt;caption&gt;Scheme of the ACM algorithm &lt;/caption&gt;<br /> &lt;/figure&gt;<br /> <br /> == Explicit/Implicit pressure calculation ==<br /> <br /> Applying divergence on \ref{NavierStokes} yields<br /> \[\nabla \cdot \frac{\partial \b{v}}{\partial t}+\nabla \cdot (\b{v}\cdot \nabla )\b{v}=-\frac{1}{\rho }{{\nabla }^{2}}p+\nabla \cdot \nu {{\nabla }^{2}}\b{v}+\nabla \cdot \b{f}\]<br /> <br /> And since $\nabla \cdot \b{v}=0$ and we can change order in $\nabla \cdot \nabla^2$ and $ \nabla^2 \cdot \nabla$ equation simplifies to<br /> \[\frac{1}{\rho }{{\nabla }^{2}}p=\nabla \cdot \b{f}-\nabla \cdot (\b{v}\cdot \nabla )\b{v}\]<br /> Now, we need boundary conditions that can be obtained by multiplying the equation with a boundary normal vector<br /> \[\b{\hat{n}}\cdot \left( \frac{\partial \b{v}}{\partial t}+(\b{v}\cdot \nabla )\b{v} \right)=\b{\hat{n}}\cdot \left( -\frac{1}{\rho }\nabla p+\nu {{\nabla }^{2}}\b{v}+\b{f} \right)\cdot \b{\hat{n}}\]<br /> \[\frac{\partial p}{\partial \b{\hat{n}}}=\left( \nu {{\nabla }^{2}}\b{v}+\b{f}-\frac{\partial \b{v}}{\partial t}-(\b{v}\cdot\nabla ) \b{v} \right)\cdot \b{\hat{n}}\]<br /> <br /> Note that using tangential boundary vector gives equivalent BCs<br /> \[\frac{\partial p}{\partial \b{\hat{t}}}=\left( \nu {{\nabla }^{2}}\b{v}+\b{f}-\frac{\partial \b{v}}{\partial t}-(\b{v}\cdot\nabla ) \b{v} \right)\cdot \b{\hat{t}}\]<br /> For no-slip boundaries BCs simplify to<br /> \[\frac{\partial p}{\partial \b{\hat{n}}}=\left( \nu {{\nabla }^{2}}\b{v}+\b{f} \right)\cdot \b{\hat{n}}\]<br /> Otherwise an appropriate expression regarding the velocity can be written, i.e. write full and taken in account velocity BCs. For example, Neumann velocity $\frac{\partial u}{\partial x}=0$ in 2D<br /> \[\frac{\partial p}{\partial x}=\left( \nu {{\nabla }^{2}}u + {{f}_{x}}-\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial y} \right)\]<br /> Note that you allready know everything about the velocity and thus you can compute all the terms explicitly.<br /> <br /> <br /> So the procedure is:<br /> <br /> * Compute Navier Stokes either explicitly or implicitly<br /> * Solve pressure equations with computed velocities<br /> * March in time<br /> <br /> Basic boundary conditions<br /> Wall: $\b{v}=0$, \[\frac{\partial p}{\partial \hat{n}}=\left( \nabla \cdot \left( \nu \nabla \b{v} \right)+\b{f} \right)\cdot \hat{n}\]<br /> Inlet: $\b{v}=\b{a}$, \[\frac{\partial p}{\partial \hat{n}}=\left( \nabla \cdot \left( \nu \nabla \b{v} \right)+\b{f}-\nabla \cdot (\rho \b{v}\b{v})-\rho \frac{\partial \b{v}}{\partial t} \right)\cdot \hat{n}\]<br /> <br /> Above system can be linearized (advection term) and solved either explicitly or implicitly.<br /> <br /> Further reading:<br /> <br /> W. D. Henshaw, A fourth-order accurate method for the incompressible Navier–Stokes equations on overlapping grids, J. Comput. Phys. 113, 13 (1994)<br /> <br /> J. C. Strikwerda, Finite difference methods for the Stokes and Navier–Stokes equations, SIAM J. Sci. Stat.<br /> Comput. 5(1), 56 (1984)<br /> <br /> == Explicit Pressure correction ==<br /> Another possibility is to solve pressure correction equation. Again Consider the momentum equation and mass continuity and discretize it explicitly<br /> \[\frac{{{\b{v}}_{2}}-{{\b{v}}_{1}}}{\Delta t}=-\frac{1}{\rho }\nabla {{p}_{1}}-({{\b{v}}_{1}}\nabla )\cdot {{\b{v}}_{1}}+\nu {{\nabla }^{2}}{{\b{v}}_{1}}+\b{f}\]<br /> Computed velocity obviously does not satisfy the mass contunity and therefore let’s call it intermediate velocity. Intermediate velocity is calculated from guessed pressure and old velocity values.<br /> \[{{\b{v}}^{inter}}=\b{v}_1 + \Delta t\left( -\frac{1}{\rho }\nabla {{p}_{1}}-({{\b{v}}_{1}}\nabla )\cdot {{\b{v}}_{1}}+\nu {{\nabla }^{2}}{{\b{v}}_{1}}+\b{f} \right)\]<br /> A correction term is added that drives velocity to divergence free field<br /> \[\nabla \cdot ({{\b{v}}^{inter}}+{{\b{v}}^{corr}})=0 \qquad \to \qquad \nabla \cdot {{\b{v}}^{inter}}=-\nabla \cdot {{\b{v}}^{corr}}\]<br /> <br /> Velocity correction is affected only by effect of pressure correction. This fact is obvious due to all terms except gradient of pressure on the right side of equation are constant.<br /> \[{{\b{v}}^{corr}}=-\frac{\Delta t}{\rho }\nabla {{p}^{corr}} \]<br /> <br /> Note that corrected velocity also satisfies boundary conditions<br /> \[\b{v}^{iter}+\b{v}^{corr}=\b{v}^{BC}\]<br /> Applying divergence and we get '''pressure correction poisson equation'''.<br /> \[\,{{\nabla }^{2}}{{p}^{corr}}\,=\frac{\rho }{\Delta t}\nabla \cdot {{\mathbf{v}}^{iter}}\,\]<br /> <br /> Boundary conditions can be obtained by mulitplying the equation with a unit normal vector $\b{\hat{n}}$<br /> \[\frac{\Delta t}{\rho }\frac{\partial {p}^{corr}}{\partial \b{\hat{n}}} = \b{\hat{n}} \cdot \left(\b{v}^{iter} - \b{v}^{BC} \right) \]<br /> The most straightforward approach, for dirichlet BCs, is to take into account velocity boundary condition in computation of intermediate velocity, and clearly in such cases, pressure boundary condition simplifies to <br /> \[\frac{\partial p^{corr}}{\partial \b{\hat{n}}} = 0 \]<br /> As ${{\b{v}}^{\operatorname{int}er}}={{\b{v}}^{BC}}$ . Another option is to explicitely compute intermediate velocity also on boundaries and then correct it through pressure correction.<br /> <br /> The pressure poisson equation is, at given boundary conditions, defined only up to a constant. One solution is to select a node and set it to a constant, e.g. p(0, 0) = 0, however much more stable approach is to enforce solution with additional condition, also referred to as a regularization<br /> \[\int_{\Omega }^{{}}{pd}\Omega =0\]<br /> \[\,{{\nabla }^{2}}{{p}^{corr}}\,-\alpha =\frac{\rho }{\Delta t}\nabla \cdot {{\b{v}}^{iter}}\,\]<br /> Where $\alpha $ stands for Lagrange multiplier. Or in discrete form<br /> \[\sum\limits_{i}{p\left( {{x}_{i}} \right)=0}\]<br /> \[\b{Mp}-\alpha \b{1}=\frac{\rho }{\Delta t}\nabla \cdot {{\b{v}}^{iter}}\]<br /> <br /> where $\b{M}$ holds Laplace shape functions, i.e. the discrete version of Laplace differential operator. <br /> <br /> Solution of a system<br /> <br /> \[\left[ \begin{matrix}<br /> {{M}_{11}} &amp; .. &amp; {{M}_{1n}} &amp; 1 \\<br /> .. &amp; .. &amp; .. &amp; 1 \\<br /> {{M}_{n1}} &amp; ... &amp; {{M}_{nn}} &amp; 1 \\<br /> 1 &amp; 1 &amp; 1 &amp; 0 \\<br /> \end{matrix} \right]\left[ \begin{matrix}<br /> {{p}_{1}} \\<br /> ... \\<br /> {{p}_{n}} \\<br /> \alpha \\<br /> \end{matrix} \right]=\frac{\rho }{\Delta t}\left[ \begin{matrix}<br /> \nabla \cdot \b{v}_{_{1}}^{\text{iter}} \\<br /> ... \\<br /> \nabla \cdot \b{v}_{n}^{\text{iter}} \\<br /> 0 \\<br /> \end{matrix} \right]\]<br /> Gives us a solution of pressure correction.<br /> <br /> == CBS Algorithm ==<br /> With explicit temporal discretization problem is formulated as<br /> \[\b{\hat{v}}={{\b{v}}_{0}}+\Delta t\left( -\nabla {{p}_{0}}+\frac{1}{Re}{{\nabla }^{2}}{{\b{v}}_{0}}-\nabla \cdot ({{\b{v}}_{0}}{{\b{v}}_{0}}) \right)\]<br /> \[p={{p}_{0}}-\xi \Delta {{t}_{F}}\nabla \b{\hat{v}}+\xi \Delta {{t}_{F}}\Delta t{{\nabla }^{2}}{{\overset{\scriptscriptstyle\frown}{P}}_{0}},\]<br /> where $\b{\hat{v}}$, $\Delta t$, $\xi$ and $\Delta t_F$ stand for intermediate velocity, time step, relaxation parameter, and artificial time step, respectively, and index 0 stands for previous time / iteration step. First, the intermediate velocity is computed from previous time step. Second, the velocity is driven towards solenoidal field by correcting the pressure. Note that no special boundary conditions for pressure are used, i.e., the pressure on boundaries is computed with the same approach as in the interior of the domain. In general, the internal iteration with an artificial time step is required until the divergence of the velocity field is not below required criteria. However, if one is interested only in a steady-state solution, the internal iteration can be skipped and $\Delta t$ equals $\Delta {{t}_{F}}$. Without internal stepping the transient of the solution is distorted by artificial compressibility effect. This approach is also known as ACM with Characteristics-based discretization of continuity equation, where the relaxation parameter relates to the artificial speed of sound [35].<br /> <br /> The relaxation parameter should be set between 1-10, lower number more stable solution.<br /> <br /> And also dimensional form<br /> <br /> \[p={{p}_{0}}-{{C}^{2}}\Delta {{t}_{F}}\rho \nabla \b{\hat{v}}+{{C}^{2}}\Delta {{t}_{F}}\Delta t{{\nabla }^{2}}{{p}_{0}},\]<br /> <br /> Where C is speed of sound [m/s]<br /> <br /> == The Stream Function - Vorticity Approach ==<br /> <br /> In two dimensions, the Navier–Stokes equations can be expressed using the<br /> stream function $\psi$ and the vorticity $\omega$ in place of the primitive variables $u$, $v$, and<br /> $p$. This involves the elimination of the pressure $p$, thus yielding one dependent<br /> variable less. In three dimensions, however, this formulation leads to six unknowns<br /> rather than four (in primitive variables), which makes this approach less attractive<br /> for that case.<br /> <br /> In the following, we briefly derive the resulting two-dimensional equations for<br /> $\psi$ and $\omega$. We begin by considering the momentum equations:<br /> \begin{align}<br /> \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}&amp; = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) + g_x \\<br /> \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}&amp; = -\frac{1}{\rho}\frac{\partial p}{\partial y} + \nu \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) + g_y<br /> \end{align}<br /> from which we can eliminate the pressure by differentiating the first equation by $y$ and the second by $x$, subtracting the second from the first and then substituting the definition for the vorticity $\omega = \partial u / \partial y - \partial v / \partial x$. In this manner we obtain the equation:<br /> \[\frac{\partial \omega}{\partial t} + \frac{\partial u}{\partial x} \omega + u \frac{\partial \omega}{\partial x} + \frac{\partial v}{\partial y}\omega + v \frac{\partial \omega}{\partial y} = \nu \left(\frac{\partial^2 \omega}{\partial x^2} + \frac{\partial^2 \omega}{\partial y^2}\right) + \left(\frac{\partial g_x}{\partial y} - \frac{\partial g_y}{\partial x}\right)\]<br /> <br /> We can get rid of the terms containing $\partial u / \partial x$ and $\partial v / \partial y$ by using the continuity equation \eqref{contuinity}. Finally using the definition for the stream function $\partial \psi / \partial y = u$ and $\partial \psi / \partial x = -v$ we can transform the above equation to what is known as the ''vorticity transport equation'':<br /> \begin{equation}<br /> \frac{\partial \omega}{\partial t} + \frac{\partial \psi}{\partial y}\frac{\partial \omega}{\partial x} - \frac{\partial \psi}{\partial x}\frac{\partial \omega}{\partial y} = \nu \Delta \omega + \left(\frac{\partial g_x}{\partial y} - \frac{\partial g_y}{\partial x}\right)<br /> \label{vorticity_transport}<br /> \end{equation}<br /> <br /> Another equation is obtained by inserting the definition of the stream function into that of the vorticity:<br /> \[\omega = \frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} = \frac{\partial}{\partial y}\left(\frac{\partial \psi}{\partial y}\right) + \frac{\partial }{\partial x}\left(\frac{\partial \psi}{\partial x}\right) = \Delta \psi\]<br /> leading to the ''Poisson equation for the stream function''<br /> \begin{equation}<br /> \Delta \psi = \omega<br /> \label{poisson_stream}<br /> \end{equation}<br /> <br /> The two equations (\ref{vorticity_transport}) and (\ref{poisson_stream}) form a nonlinear coupled system of equations in which pressure has been eliminated and in which the continuity equation has been satisfied automatically.<br /> <br /> = Numerical examples =<br /> * [[Lid driven cavity]]<br /> * [[Burger's equation]]<br /> * [[de Vahl Davis natural convection test]]<br /> * [[Natural convection in 3D irregular domain]]<br /> * [[Natural convection from heated cylinder]]<br /> * [[Natural convection between concentric cylinders]]<br /> * [[Non-Newtonian fluid]]</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Burgers%27_equation&diff=3560 2023-12-27T12:31:16Z

<p>Zigavaupotic: Zigavaupotic moved page Burgers equation to Burger's equation</p> <hr /> <div>Burgers equation</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Burgers_equation&diff=3561 2023-12-27T12:31:16Z

<p>Zigavaupotic: Zigavaupotic moved page Burgers equation to Burger's equation</p> <hr /> <div>#REDIRECT [[Burger's equation]]</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Fluid_Mechanics&diff=3559 2023-12-27T12:30:28Z

<p>Zigavaupotic: </p> <hr /> <div>= Introduction =<br /> Computational fluid dynamics (CFD) is a field of a great interest among researchers in many fields of science, e.g. studying mathematical fundaments of numerical methods, developing novel physical models, improving computer implementations, and many others. Pushing the limits of all involved fields of science helps community to deepen the understanding of several natural and technological phenomena. Weather forecast, ocean dynamics, water transport, casting, various energetic studies, etc., are just few examples where fluid dynamics plays a crucial role. The core problem of the CFD is solving the Navier-Stokes Equation or its variants, e.g. Darcy or Brinkman equation for flow in porous media. Here, we discuss basic algorithms for solving CFD problems. Check reference list on the [[Medusa|main page]] for more details about related work.<br /> <br /> Long story short, we want to solve<br /> \begin{equation}<br /> \frac{\partial \b{v}}{\partial t}+\nabla \cdot ( \rho \b{vv}) = -\nabla p+ \nabla(\mu \nabla\b{v})+\b{f}<br /> \label{NavierStokes}<br /> \end{equation}<br /> also known as a Navier-Stokes equation. <br /> <br /> Note that the $\b{v}\b{v}$ stands for the tensor or dyadic product \[ \b{v}\b{v} = \b{v}\otimes\b{v} = \b{v}\b{v}^\T = \left[ \begin{matrix}<br /> {{v}_{1}}{{v}_{1}} &amp; \cdots &amp; {{v}_{1}}{{v}_{n}} \\<br /> \vdots &amp; \ddots &amp; \vdots \\<br /> {{v}_{n}}{{v}_{1}} &amp; \cdots &amp; {{v}_{n}}{{v}_{n}} \\<br /> \end{matrix} \right]\]<br /> <br /> Combining equation \ref{NavierStokes} and mass contiunity equation<br /> \begin{equation}<br /> \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\bf{u}}) = 0.<br /> \end{equation}<br /> <br /> simplifies advection term into (see [https://en.wikipedia.org/wiki/Derivation_of_the_Navier%E2%80%93Stokes_equations] for derivation)<br /> <br /> \[\frac{\partial \left( \rho \b{v} \right)}{\partial t}+\nabla \cdot \left( \rho \b{vv} \right)=\frac{\partial \left( \rho \b{v} \right)}{\partial t}+(\rho \b{v}\cdot \nabla )\cdot \b{v}. \]<br /> <br /> An example of advection term in 2D would therefore be<br /> \[\left( \b{v}\cdot \nabla \right)\b{v}=\left( \left( \begin{matrix}<br /> u \\<br /> v \\<br /> \end{matrix} \right) \cdot \left( \begin{matrix}<br /> \frac{\partial }{\partial x} \\<br /> \frac{\partial }{\partial y} \\<br /> \end{matrix} \right) \right)\left( \begin{matrix}<br /> u \\<br /> v \\<br /> \end{matrix} \right)=\left( \begin{matrix}<br /> u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} \\<br /> u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} \\<br /> \end{matrix} \right)\]<br /> <br /> <br /> In many cases we are interested in the incompressible fluids (Ma&lt;0.3), reducing the continuity equation to<br /> \begin{equation}<br /> \nabla \cdot \b{v}=0<br /> \label{contuinity}<br /> \end{equation}<br /> <br /> <br /> The goal of CFD is to solve system \ref{NavierStokes} and \ref{contuinity}. It is obvious that a special treatment will be needed to couple both equations. In the following discussion we cover some basic approaches, how this can be accomplished.<br /> <br /> = Solutions algorithms =<br /> == Artificial compressibility method ==<br /> The simplest, completely explicit approach, is an artificial compressibility method (ACM), where a compressibility term is included in the mass continuity<br /> \[\frac{\partial \b{v}}{\partial t}+(\b{v}\cdot\nabla )\b{v}=-\frac{1}{\rho }\nabla p+\nu {{\nabla }^{2}}\b{v}+\b{f}\]<br /> \[\frac{ 1 }{ \rho } \frac{\partial \rho }{\partial t}+\nabla \cdot \b{v}=0\]<br /> \[\frac{ 1 }{ \rho } \frac{\partial \rho }{\partial p}\frac{\partial p}{\partial t}+\nabla \cdot \b{v}=0\]<br /> Now, the above system can be solved directly.<br /> <br /> The addition of the time derivative of the pressure term physically means that waves of finite speed (the propagation of which depends on the magnitude of the ACM)<br /> are introduced into the flow field as a mean to distribute the pressure within the domain. In a true<br /> incompressible flow, the pressure field is affected instantaneously throughout the whole domain. In ACM there is a time delay between the flow disturbance and its effect on the<br /> pressure field. Upon rearranging the equation yields<br /> \[\frac{\partial p}{\partial t}+\rho {{C}^{2}}\nabla \cdot \b{v}=0\]<br /> where the continuity equation is perturbed by the quantity $\frac{\partial p}{\partial t}$ denominated herein<br /> as the AC parameter/artificial sound speed recognized by<br /> $C$ [m/s] - speed of sound<br /> \[\frac{1}{C^2}=\frac{\partial \rho }{\partial p}\]<br /> Or in other words<br /> \[C^2=\left( \frac{\partial p}{\partial \rho}\right)_S\]<br /> where $\rho$ is the density of the material. It follows, by replacing partial derivatives, that the isentropic compressibility can be expressed as:<br /> \[\beta =\frac{1}{\rho {{C}^{2}}}\]<br /> The evaluation of the local ACM parameter in incompressible flows is inspired by the<br /> speed of sound computations in compressible flows (for instance, from the perfect gas law).<br /> However, in the incompressible flow situation, employing such a relation is difficult, but an artificial<br /> relation can be developed from the convective and diffusive velocities.<br /> Reverting to the justification of continuity modification, it can be immediately seen that the<br /> artificial sound speed must be sufficiently large to have a significant regularizing effect and at<br /> the same time must be as small as possible to minimizing perturbations on the incompressibility<br /> equation. Therefore, $C$ influences the convergence rate and stability of the solution method. In other words,<br /> assists in reducing large disparity in the eigenvalues, leading to a well-conditioned system. <br /> The $C$ can be '''estimated''' with<br /> <br /> \[ C = \beta \max \left( \left|\b{v}\right|_2, \left|\b{v}_{ref}\right|_2 \right),\]<br /> where $\b{v}_{ref}$ stands for a reference velocity. <br /> Values for $\beta$ in the range of 1–10 are recommended for better convergence to the steady state at which the<br /> mass conservation is enforced. In addition, Equation ensures that $C$ does not reach zero at stagnation points<br /> that cause instabilities in pseudo-time, effecting convergence<br /> <br /> Note, that for more complex simulation an internal iteration loops is required before marching in time. <br /> <br /> <br /> &lt;figure id=&quot;ACM&quot;&gt;<br /> [[File:ACM BlockDiagram.png|600px]]<br /> &lt;caption&gt;Scheme of the ACM algorithm &lt;/caption&gt;<br /> &lt;/figure&gt;<br /> <br /> == Explicit/Implicit pressure calculation ==<br /> <br /> Applying divergence on \ref{NavierStokes} yields<br /> \[\nabla \cdot \frac{\partial \b{v}}{\partial t}+\nabla \cdot (\b{v}\cdot \nabla )\b{v}=-\frac{1}{\rho }{{\nabla }^{2}}p+\nabla \cdot \nu {{\nabla }^{2}}\b{v}+\nabla \cdot \b{f}\]<br /> <br /> And since $\nabla \cdot \b{v}=0$ and we can change order in $\nabla \cdot \nabla^2$ and $ \nabla^2 \cdot \nabla$ equation simplifies to<br /> \[\frac{1}{\rho }{{\nabla }^{2}}p=\nabla \cdot \b{f}-\nabla \cdot (\b{v}\cdot \nabla )\b{v}\]<br /> Now, we need boundary conditions that can be obtained by multiplying the equation with a boundary normal vector<br /> \[\b{\hat{n}}\cdot \left( \frac{\partial \b{v}}{\partial t}+(\b{v}\cdot \nabla )\b{v} \right)=\b{\hat{n}}\cdot \left( -\frac{1}{\rho }\nabla p+\nu {{\nabla }^{2}}\b{v}+\b{f} \right)\cdot \b{\hat{n}}\]<br /> \[\frac{\partial p}{\partial \b{\hat{n}}}=\left( \nu {{\nabla }^{2}}\b{v}+\b{f}-\frac{\partial \b{v}}{\partial t}-(\b{v}\cdot\nabla ) \b{v} \right)\cdot \b{\hat{n}}\]<br /> <br /> Note that using tangential boundary vector gives equivalent BCs<br /> \[\frac{\partial p}{\partial \b{\hat{t}}}=\left( \nu {{\nabla }^{2}}\b{v}+\b{f}-\frac{\partial \b{v}}{\partial t}-(\b{v}\cdot\nabla ) \b{v} \right)\cdot \b{\hat{t}}\]<br /> For no-slip boundaries BCs simplify to<br /> \[\frac{\partial p}{\partial \b{\hat{n}}}=\left( \nu {{\nabla }^{2}}\b{v}+\b{f} \right)\cdot \b{\hat{n}}\]<br /> Otherwise an appropriate expression regarding the velocity can be written, i.e. write full and taken in account velocity BCs. For example, Neumann velocity $\frac{\partial u}{\partial x}=0$ in 2D<br /> \[\frac{\partial p}{\partial x}=\left( \nu {{\nabla }^{2}}u + {{f}_{x}}-\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial y} \right)\]<br /> Note that you allready know everything about the velocity and thus you can compute all the terms explicitly.<br /> <br /> <br /> So the procedure is:<br /> <br /> * Compute Navier Stokes either explicitly or implicitly<br /> * Solve pressure equations with computed velocities<br /> * March in time<br /> <br /> Basic boundary conditions<br /> Wall: $\b{v}=0$, \[\frac{\partial p}{\partial \hat{n}}=\left( \nabla \cdot \left( \nu \nabla \b{v} \right)+\b{f} \right)\cdot \hat{n}\]<br /> Inlet: $\b{v}=\b{a}$, \[\frac{\partial p}{\partial \hat{n}}=\left( \nabla \cdot \left( \nu \nabla \b{v} \right)+\b{f}-\nabla \cdot (\rho \b{v}\b{v})-\rho \frac{\partial \b{v}}{\partial t} \right)\cdot \hat{n}\]<br /> <br /> Above system can be linearized (advection term) and solved either explicitly or implicitly.<br /> <br /> Further reading:<br /> <br /> W. D. Henshaw, A fourth-order accurate method for the incompressible Navier–Stokes equations on overlapping grids, J. Comput. Phys. 113, 13 (1994)<br /> <br /> J. C. Strikwerda, Finite difference methods for the Stokes and Navier–Stokes equations, SIAM J. Sci. Stat.<br /> Comput. 5(1), 56 (1984)<br /> <br /> == Explicit Pressure correction ==<br /> Another possibility is to solve pressure correction equation. Again Consider the momentum equation and mass continuity and discretize it explicitly<br /> \[\frac{{{\b{v}}_{2}}-{{\b{v}}_{1}}}{\Delta t}=-\frac{1}{\rho }\nabla {{p}_{1}}-({{\b{v}}_{1}}\nabla )\cdot {{\b{v}}_{1}}+\nu {{\nabla }^{2}}{{\b{v}}_{1}}+\b{f}\]<br /> Computed velocity obviously does not satisfy the mass contunity and therefore let’s call it intermediate velocity. Intermediate velocity is calculated from guessed pressure and old velocity values.<br /> \[{{\b{v}}^{inter}}=\b{v}_1 + \Delta t\left( -\frac{1}{\rho }\nabla {{p}_{1}}-({{\b{v}}_{1}}\nabla )\cdot {{\b{v}}_{1}}+\nu {{\nabla }^{2}}{{\b{v}}_{1}}+\b{f} \right)\]<br /> A correction term is added that drives velocity to divergence free field<br /> \[\nabla \cdot ({{\b{v}}^{inter}}+{{\b{v}}^{corr}})=0 \qquad \to \qquad \nabla \cdot {{\b{v}}^{inter}}=-\nabla \cdot {{\b{v}}^{corr}}\]<br /> <br /> Velocity correction is affected only by effect of pressure correction. This fact is obvious due to all terms except gradient of pressure on the right side of equation are constant.<br /> \[{{\b{v}}^{corr}}=-\frac{\Delta t}{\rho }\nabla {{p}^{corr}} \]<br /> <br /> Note that corrected velocity also satisfies boundary conditions<br /> \[\b{v}^{iter}+\b{v}^{corr}=\b{v}^{BC}\]<br /> Applying divergence and we get '''pressure correction poisson equation'''.<br /> \[\,{{\nabla }^{2}}{{p}^{corr}}\,=\frac{\rho }{\Delta t}\nabla \cdot {{\mathbf{v}}^{iter}}\,\]<br /> <br /> Boundary conditions can be obtained by mulitplying the equation with a unit normal vector $\b{\hat{n}}$<br /> \[\frac{\Delta t}{\rho }\frac{\partial {p}^{corr}}{\partial \b{\hat{n}}} = \b{\hat{n}} \cdot \left(\b{v}^{iter} - \b{v}^{BC} \right) \]<br /> The most straightforward approach, for dirichlet BCs, is to take into account velocity boundary condition in computation of intermediate velocity, and clearly in such cases, pressure boundary condition simplifies to <br /> \[\frac{\partial p^{corr}}{\partial \b{\hat{n}}} = 0 \]<br /> As ${{\b{v}}^{\operatorname{int}er}}={{\b{v}}^{BC}}$ . Another option is to explicitely compute intermediate velocity also on boundaries and then correct it through pressure correction.<br /> <br /> The pressure poisson equation is, at given boundary conditions, defined only up to a constant. One solution is to select a node and set it to a constant, e.g. p(0, 0) = 0, however much more stable approach is to enforce solution with additional condition, also referred to as a regularization<br /> \[\int_{\Omega }^{{}}{pd}\Omega =0\]<br /> \[\,{{\nabla }^{2}}{{p}^{corr}}\,-\alpha =\frac{\rho }{\Delta t}\nabla \cdot {{\b{v}}^{iter}}\,\]<br /> Where $\alpha $ stands for Lagrange multiplier. Or in discrete form<br /> \[\sum\limits_{i}{p\left( {{x}_{i}} \right)=0}\]<br /> \[\b{Mp}-\alpha \b{1}=\frac{\rho }{\Delta t}\nabla \cdot {{\b{v}}^{iter}}\]<br /> <br /> where $\b{M}$ holds Laplace shape functions, i.e. the discrete version of Laplace differential operator. <br /> <br /> Solution of a system<br /> <br /> \[\left[ \begin{matrix}<br /> {{M}_{11}} &amp; .. &amp; {{M}_{1n}} &amp; 1 \\<br /> .. &amp; .. &amp; .. &amp; 1 \\<br /> {{M}_{n1}} &amp; ... &amp; {{M}_{nn}} &amp; 1 \\<br /> 1 &amp; 1 &amp; 1 &amp; 0 \\<br /> \end{matrix} \right]\left[ \begin{matrix}<br /> {{p}_{1}} \\<br /> ... \\<br /> {{p}_{n}} \\<br /> \alpha \\<br /> \end{matrix} \right]=\frac{\rho }{\Delta t}\left[ \begin{matrix}<br /> \nabla \cdot \b{v}_{_{1}}^{\text{iter}} \\<br /> ... \\<br /> \nabla \cdot \b{v}_{n}^{\text{iter}} \\<br /> 0 \\<br /> \end{matrix} \right]\]<br /> Gives us a solution of pressure correction.<br /> <br /> == CBS Algorithm ==<br /> With explicit temporal discretization problem is formulated as<br /> \[\b{\hat{v}}={{\b{v}}_{0}}+\Delta t\left( -\nabla {{p}_{0}}+\frac{1}{Re}{{\nabla }^{2}}{{\b{v}}_{0}}-\nabla \cdot ({{\b{v}}_{0}}{{\b{v}}_{0}}) \right)\]<br /> \[p={{p}_{0}}-\xi \Delta {{t}_{F}}\nabla \b{\hat{v}}+\xi \Delta {{t}_{F}}\Delta t{{\nabla }^{2}}{{\overset{\scriptscriptstyle\frown}{P}}_{0}},\]<br /> where $\b{\hat{v}}$, $\Delta t$, $\xi$ and $\Delta t_F$ stand for intermediate velocity, time step, relaxation parameter, and artificial time step, respectively, and index 0 stands for previous time / iteration step. First, the intermediate velocity is computed from previous time step. Second, the velocity is driven towards solenoidal field by correcting the pressure. Note that no special boundary conditions for pressure are used, i.e., the pressure on boundaries is computed with the same approach as in the interior of the domain. In general, the internal iteration with an artificial time step is required until the divergence of the velocity field is not below required criteria. However, if one is interested only in a steady-state solution, the internal iteration can be skipped and $\Delta t$ equals $\Delta {{t}_{F}}$. Without internal stepping the transient of the solution is distorted by artificial compressibility effect. This approach is also known as ACM with Characteristics-based discretization of continuity equation, where the relaxation parameter relates to the artificial speed of sound [35].<br /> <br /> The relaxation parameter should be set between 1-10, lower number more stable solution.<br /> <br /> And also dimensional form<br /> <br /> \[p={{p}_{0}}-{{C}^{2}}\Delta {{t}_{F}}\rho \nabla \b{\hat{v}}+{{C}^{2}}\Delta {{t}_{F}}\Delta t{{\nabla }^{2}}{{p}_{0}},\]<br /> <br /> Where C is speed of sound [m/s]<br /> <br /> == The Stream Function - Vorticity Approach ==<br /> <br /> In two dimensions, the Navier–Stokes equations can be expressed using the<br /> stream function $\psi$ and the vorticity $\omega$ in place of the primitive variables $u$, $v$, and<br /> $p$. This involves the elimination of the pressure $p$, thus yielding one dependent<br /> variable less. In three dimensions, however, this formulation leads to six unknowns<br /> rather than four (in primitive variables), which makes this approach less attractive<br /> for that case.<br /> <br /> In the following, we briefly derive the resulting two-dimensional equations for<br /> $\psi$ and $\omega$. We begin by considering the momentum equations:<br /> \begin{align}<br /> \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}&amp; = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) + g_x \\<br /> \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}&amp; = -\frac{1}{\rho}\frac{\partial p}{\partial y} + \nu \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) + g_y<br /> \end{align}<br /> from which we can eliminate the pressure by differentiating the first equation by $y$ and the second by $x$, subtracting the second from the first and then substituting the definition for the vorticity $\omega = \partial u / \partial y - \partial v / \partial x$. In this manner we obtain the equation:<br /> \[\frac{\partial \omega}{\partial t} + \frac{\partial u}{\partial x} \omega + u \frac{\partial \omega}{\partial x} + \frac{\partial v}{\partial y}\omega + v \frac{\partial \omega}{\partial y} = \nu \left(\frac{\partial^2 \omega}{\partial x^2} + \frac{\partial^2 \omega}{\partial y^2}\right) + \left(\frac{\partial g_x}{\partial y} - \frac{\partial g_y}{\partial x}\right)\]<br /> <br /> We can get rid of the terms containing $\partial u / \partial x$ and $\partial v / \partial y$ by using the continuity equation \eqref{contuinity}. Finally using the definition for the stream function $\partial \psi / \partial y = u$ and $\partial \psi / \partial x = -v$ we can transform the above equation to what is known as the ''vorticity transport equation'':<br /> \begin{equation}<br /> \frac{\partial \omega}{\partial t} + \frac{\partial \psi}{\partial y}\frac{\partial \omega}{\partial x} - \frac{\partial \psi}{\partial x}\frac{\partial \omega}{\partial y} = \nu \Delta \omega + \left(\frac{\partial g_x}{\partial y} - \frac{\partial g_y}{\partial x}\right)<br /> \label{vorticity_transport}<br /> \end{equation}<br /> <br /> Another equation is obtained by inserting the definition of the stream function into that of the vorticity:<br /> \[\omega = \frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} = \frac{\partial}{\partial y}\left(\frac{\partial \psi}{\partial y}\right) + \frac{\partial }{\partial x}\left(\frac{\partial \psi}{\partial x}\right) = \Delta \psi\]<br /> leading to the ''Poisson equation for the stream function''<br /> \begin{equation}<br /> \Delta \psi = \omega<br /> \label{poisson_stream}<br /> \end{equation}<br /> <br /> The two equations (\ref{vorticity_transport}) and (\ref{poisson_stream}) form a nonlinear coupled system of equations in which pressure has been eliminated and in which the continuity equation has been satisfied automatically.<br /> <br /> = Numerical examples =<br /> * [[Lid driven cavity]]<br /> * [[Burgers equation]]<br /> * [[de Vahl Davis natural convection test]]<br /> * [[Natural convection in 3D irregular domain]]<br /> * [[Natural convection from heated cylinder]]<br /> * [[Natural convection between concentric cylinders]]<br /> * [[Non-Newtonian fluid]]</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Burgers%27_equation&diff=3558 2023-12-27T12:27:49Z

<p>Zigavaupotic: Created page with &quot;Burgers equation&quot;</p> <hr /> <div>Burgers equation</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Solving_system&diff=3557 2023-12-27T12:07:21Z

<p>Zigavaupotic: /* ordinary least squares */</p> <hr /> <div>= Basics of solving square linear systems =<br /> <br /> Solving a system $Ax = b$, where $A$ is a $n\times n$ matrix of scalars (complex, real, or otherwise) is possible if $A$ is invertible, and the solution is given theoretically by $x = A^{-1}b$.<br /> <br /> <br /> Let $\hat{x}$ be the solution of $\hat{A} \hat{x} = \hat{b}$, where the matrix $A$ is disturbed as $\hat{A} = A + \Delta{A}$ and the right hand side as $\hat{b} = b + \Delta b$.<br /> <br /> The absolute error of $x$ is $\Delta x = \|\hat{x} - x\|$ and the relative errors are defined as $\delta A = \|\Delta A\|/\|A\|$ and $\delta b = \|\Delta b\|/\|b\|$ and $\delta x = \|\Delta x\|/\|x\|$.<br /> We assume that relative errors in inputs are reasonably small.<br /> <br /> <br /> If only $b$ is disturbed, the error is given by <br /> $$\delta x \leq \kappa(A) \delta b,$$ where $\kappa(A) = \|A\|\|A^{-1}\|$ is the ''condition number'' of matrix $A$.<br /> <br /> In general, the error is given by<br /> <br /> $$\delta x \leq \frac{\kappa(A)}{1 - \kappa(A) \delta A} (\delta A + \delta b).$$<br /> <br /> Therefore, if the condition number $\kappa(A)$ is large, the accuracy of the solution of the system can deteriorate.<br /> <br /> If $A$ is not invertible, Moore-Penrose pseudoinverse can be used to compute a &quot;solution&quot;: $x = A^+b$.<br /> <br /> = Basics of solving underdetermined linear systems = <br /> = Least norm solution =<br /> Consider the system $Ax = b$, where $A$ is a $n \times m$ matrix of full rank and $n \leq m$. The number of variables $m$ is larger than the number of equations. This means that all equations can be satisfied <br /> and some degrees of freedom are left to obtain various properties of the solution. The common problem is to find the least-norm solution: find $x$, such that<br /> $$\min_x \|x\|_2, \text{ such that } Ax = b.$$<br /> <br /> The solution is given by $x = A^\T (AA^\T)^{-1} b$, and it can be obtained by solving the $n \times n$ system $AA^\T y = b$ and computing $x = A^\T y$. <br /> The solution can also be written as $x = A^+ b$. <br /> <br /> To validate the correctness, consider any other $z$, such that $Az=b$. Then, $$<br /> (z-x)^\T x = (z-x)^\T A^\T (AA^\T)^{-1} b = (A(x-z))^\T (AA^\T)^{-1} b = (b-b)^\T (AA^\T)^{-1} b = 0,<br /> $$<br /> and consequently<br /> $$<br /> \|z\|^2 = \|z-x+x\|^2 = \|z-x\|^2 + 2(z-x)^\T x + \|x\|^2 = \|z-x\|^2 + \|x\|^2 \geq \|x\|^2,<br /> $$<br /> which shows that $x$ has smaller norm than any other solution $z$.<br /> <br /> Alternative derivation is to minimize $x^\T x$ with respect to $Ax=b$ using Lagrangian multipliers.<br /> <br /> = Weighted case =<br /> <br /> Weighted case is similar, except that the minimization problem becomes $$\min_x \|x\|_{2,w}, \text{ such that } Ax = b.$$<br /> The minimization can be written $\min_x \|Wx\|_{2}, \text{ such that } Ax = b$, or substituting $Wx = y$, we obtain the following<br /> ordinary least norm problem:<br /> $$\min_y \|y\|_{2}, \text{ such that } AW^{-1}y = b,$$<br /> with matrix $AW^{-1}$ and right hand side $b$, and with solution $x = W^{-1} y$.<br /> <br /> Thus, the final solution is given as $x = W^{-2} A^\T (A W^{-2} A^\T)^{-1} b$ and can be computed by first solving <br /> $n \times n$ system $AW^{-2}A^\T y = b$ and computing $x = W^{-2} A^\T y$.<br /> <br /> = Numerical solving =<br /> <br /> Techniques for numerical solving of (weighted) underdetermined systems are described below. As list of decompositions available in Eigen is available here:<br /> https://eigen.tuxfamily.org/dox/group__LeastSquares.html<br /> <br /> As when solving overdetermined systems if a lot of systems $Ax = b_i$ need to be solved, which only differ in the right hand side, the matrix $A$ can only be decomposed once, with<br /> subsequent solutions taking $O(n^2)$ time.<br /> <br /> == Normal equations ==<br /> <br /> We can solve the system directly, by solving the $n \times n$ system $AA^\T y = b$ using Cholesky, and computing $x = A^\T y$.<br /> <br /> This decomposition takes $\frac13 n^3 + 2n^2m$ time.<br /> <br /> == QR decomposition ==<br /> <br /> Writing $A^\T = Q_1R_1$ as when solving overdetermined systems, we can get the solution<br /> $$x = Q_1R_1 (R_1^\T Q_1^\T Q_1 R_1)^{-1}b = Q_1 R_1 R_1^{-1} R_1^{-\T} b = Q_1 R_1^{-\T} b.$$<br /> <br /> The computational costs are the same as in the overdetermined case.<br /> <br /> == SVD decomposition ==<br /> <br /> The pseudoinverse can be computed using SVD and the solution obtained as $x = A^{+} b$.<br /> <br /> The computational costs are the same as in the overdetermined case.<br /> <br /> == Decompositions ==<br /> <br /> Solving linear systems is not done by computing inverses of $A$ but by suitably decomposing $A$. <br /> Many decompositions are available, with $LU$ decomposition with partial pivoting being the go-to option.<br /> with $\frac23 n^3$ cost. Additionally, we have $QR$ and $SVD$ decompositions if more robust decompositions are desired, at a greater cost.<br /> <br /> See full the list of Eigen's available decompositions: https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html<br /> <br /> If many systems $Ax = b_i$ need to be saved with only different right hand sides, the decomposition of $A$ can be computed in $O(n^3)$ and<br /> stored with subsequent solutions only taking $O(n^2)$ time.<br /> <br /> If the matrix is symmetric and positive definite, Cholesky (called $VV^\T$ or $LL^\T$) decomposition can be used, with $\frac13 n^3$ cost.<br /> <br /> = Basics of solving overdetermined linear systems =<br /> == Ordinary least squares ==<br /> <br /> Let $A$ be a matrix of size $n \times m$ and $b$ be a vector of size $n$, with $n \geq m$ and let $A$ have full rank. If $n=m$ the system $Ax=b$ is square, and can be treated as such, however the<br /> derivation below works as well.<br /> <br /> The system $Ax=b$ has $n$ equations for $m$ unknowns $x_1, \ldots, x_m$ and is overdetermined. The solution is most often sought in the least-squares sense, defining the solution $x$ to be the solution of the following minimisation problem<br /> $$ \min_x \|Ax-b\|_2. $$<br /> The solution can be obtained using calculus, by looking for stationary points of function $\|Ax-b\|_2^2 = (Ax-b)^\T(Ax-b)$ (which has the same minimum), and getting the system of<br /> normal equations $$A^\T A x = A^\T b.$$ The solution can be expressed theoretically as<br /> $x = (A^\T A)^{-1} A^\T b$.<br /> <br /> == Weighted least squares ==<br /> <br /> The weighted least squares problem has the same setup as the ordinary least squares, but we instead minimise the norm $\|Ax-b\|_{2,w}$, where $\|z\|_{2,w} = \sum_{i=1}^n (w_iz_i)^2$. We assume all weights $w_i$ are positive.<br /> <br /> By observing that $\|z\|_{2,w} = \|Wz\|_2$, where $W = \operatorname{diag}(w_1, \ldots, w_n)$, we can reduce the problem to ordinary least squares.<br /> Writing $\|Ax-b\|_{2,w} = \|W(Ax-b)\|_2 = \|WA x - Wb\|_2$, we obtain an ordinary least squares problem with matrix $WA$ and vector $Wb$, which has the <br /> normal equations $$A^\T W^2Ax = A^\T W^2b$$ and the solution $x = (A^\T W^2A)^{-1} A^\T W^2b$.<br /> <br /> = Numerical solving =<br /> <br /> Techniques for numerical solving of (weighted) least squares systems are described below. As list of decompositions available in Eigen is available here:<br /> https://eigen.tuxfamily.org/dox/group__LeastSquares.html<br /> <br /> As when solving linear systems if a lot of systems $Ax = b_i$ need to be solved, which only differ in the right hand side, the matrix $A$ can only be decomposed once, with<br /> subsequent solutions taking $O(n^2)$ time.<br /> <br /> == Normal equations ==<br /> The normal equations above can be directly used to solve the system. The matrix $A^\T A$ is symmetric and positive definite, and<br /> normal equations can be solved using Cholesky decomposition. However, this can be <br /> poorly behaved with respect to round-off errors since the condition number $\kappa(A^\T A)$ is the square<br /> of $\kappa(A)$.<br /> <br /> Complexity of Cholesky decomposition is $\frac{m^3}{3}$ and complexity of matrix multiplication $2m^2n$. To preform the Cholesky decomposition the rank of $A$ must be full.<br /> <br /> '''Pros:'''<br /> * simple to implement<br /> * low computational complexity<br /> <br /> '''Cons:'''<br /> * numerically unstable<br /> * full rank requirement<br /> <br /> == [https://en.wikipedia.org/wiki/QR_decomposition $QR$ Decomposition] ==<br /> The square $n\times m$ matrix $A$ is decomposed as <br /> $$<br /> A = QR = <br /> \begin{bmatrix} Q_1 &amp; Q_2 \end{bmatrix} <br /> \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1R_1,$$<br /> where <br /> $Q$ is a unitary matrix ($Q^{-1}=Q^T$) of size $n\times n$, $R$ trapezoidal and is the same size as $A$ ($n \times m$),<br /> $Q_1$ is also the same as $A$ ($n \times m$) and $R_1$ is a triangular $m \times m$ matrix.<br /> A useful property of a unitary matrices is that multiplying with them does not alter the (Euclidean) norm of a vector, i.e.,<br /> $$ \| Qx \| = \| x \|. $$<br /> Therefore we can say<br /> $$<br /> \| Ax- b \| = \| Q^\T (Ax-b) \| = \| Q^\T Ax - Q^\T b \|<br /> = \| Q^\T QRx - Q^\T b \| = \| [R_1; 0] x - [Q_1, Q_2]^\T b \|<br /> = \| R_1x - Q_1^\T b \| + \| Q_2^\T b \|<br /> $$<br /> Of the two terms, we have no control over the second, and we can render the first one<br /> zero by solving<br /> $$ R_1x = Q_1^\T b, $$<br /> which results in a minimum $x = R_1^{-1} Q_1^\T b$.<br /> <br /> When computing the $QR$ decomposition, only the ''thin'' decomposition $Q_1$, $R_1$ is needed to solve the least squares problem.<br /> <br /> The complexity of $QR$ decomposition is $2nm^2−2m^3/3$.<br /> <br /> &lt;strong&gt;Pros:&lt;/strong&gt; better stability in comparison to normal equations, &lt;strong&gt; cons: &lt;/strong&gt;higher complexity<br /> <br /> == [https://en.wikipedia.org/wiki/Singular_value_decomposition SVD decomposition] ==<br /> In linear algebra, the [https://en.wikipedia.org/wiki/Singular_value_decomposition singular value decomposition (SVD)]<br /> is a factorization of a real or complex matrix. <br /> <br /> The singular value decomposition of an $n \times m$ real or complex<br /> matrix $A$ is a factorization of the form $A= U\Sigma V^\T$, where<br /> $U$ is an $n \times n$ unitary matrix, $\Sigma$ is an $n \times m$<br /> rectangular diagonal matrix with non-negative real numbers on the diagonal, and<br /> $V$ is an $m \times m$ unitary matrix. The diagonal entries<br /> $\Sigma_{ii}$ are known as the singular values of $A$. The $n$ columns of<br /> $U$ and the $n$ columns of $V$ are called the left-singular vectors and<br /> right-singular vectors of $A$, respectively.<br /> <br /> The singular value decomposition and the eigendecomposition are closely<br /> related. Namely:<br /> <br /> * The left-singular vectors of $A$ are eigen vectors of $A A^\T$.<br /> * The right-singular vectors of $A$ are eigen vectors of $A^\T A$.<br /> * The non-zero singular values of $A$ (found on the diagonal entries of $\Sigma$) are the square roots of the non-zero eigenvalues of both $A^\T A$ and $A^\T A$.<br /> <br /> with SVD we can write $A$ as \[A= U\Sigma V^\T \] where are $U$ and $V$ again unitary matrices and $\Sigma$<br /> stands for diagonal matrix of singular values.<br /> <br /> SVD can be used to compute the pseudoinverse<br /> \[ A^{+} = \left( U\Sigma V^\T\right)^{-1} = V \Sigma^{-1}U^\T \]<br /> where $\Sigma^{-1}$ is trivial, just replace every non-zero diagonal entry by<br /> its reciprocal and transpose the resulting matrix. <br /> <br /> This can be used to solve the least squares problem, ordinary systems or underdetermined systems.<br /> Note that usually only the ''thin'' variant is needed.<br /> <br /> One can also choose a threshold $t$<br /> below which the singular value is considered as $0$ and truncate all <br /> singular values, whose magnitude is below $t$. Because we do not invert very small values,<br /> stability of the solution increases, at the cost of a slightly inexact solution.<br /> This approach is called &quot;truncated SVD decomposition&quot;, which is a form of regularization of badly conditioned problems.<br /> Another approach would be [https://en.wikipedia.org/wiki/Tikhonov_regularization Tikhonov regularization].<br /> <br /> SVD decomposition complexity \[ 2nm^2+2m^3 = O(n^3) \]<br /> <br /> &lt;strong&gt;Pros:&lt;/strong&gt; stable, &lt;strong&gt; cons: &lt;/strong&gt;high complexity</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Solving_system&diff=3556 2023-12-27T12:06:59Z

<p>Zigavaupotic: </p> <hr /> <div>= Basics of solving square linear systems =<br /> <br /> Solving a system $Ax = b$, where $A$ is a $n\times n$ matrix of scalars (complex, real, or otherwise) is possible if $A$ is invertible, and the solution is given theoretically by $x = A^{-1}b$.<br /> <br /> <br /> Let $\hat{x}$ be the solution of $\hat{A} \hat{x} = \hat{b}$, where the matrix $A$ is disturbed as $\hat{A} = A + \Delta{A}$ and the right hand side as $\hat{b} = b + \Delta b$.<br /> <br /> The absolute error of $x$ is $\Delta x = \|\hat{x} - x\|$ and the relative errors are defined as $\delta A = \|\Delta A\|/\|A\|$ and $\delta b = \|\Delta b\|/\|b\|$ and $\delta x = \|\Delta x\|/\|x\|$.<br /> We assume that relative errors in inputs are reasonably small.<br /> <br /> <br /> If only $b$ is disturbed, the error is given by <br /> $$\delta x \leq \kappa(A) \delta b,$$ where $\kappa(A) = \|A\|\|A^{-1}\|$ is the ''condition number'' of matrix $A$.<br /> <br /> In general, the error is given by<br /> <br /> $$\delta x \leq \frac{\kappa(A)}{1 - \kappa(A) \delta A} (\delta A + \delta b).$$<br /> <br /> Therefore, if the condition number $\kappa(A)$ is large, the accuracy of the solution of the system can deteriorate.<br /> <br /> If $A$ is not invertible, Moore-Penrose pseudoinverse can be used to compute a &quot;solution&quot;: $x = A^+b$.<br /> <br /> = Basics of solving underdetermined linear systems = <br /> = Least norm solution =<br /> Consider the system $Ax = b$, where $A$ is a $n \times m$ matrix of full rank and $n \leq m$. The number of variables $m$ is larger than the number of equations. This means that all equations can be satisfied <br /> and some degrees of freedom are left to obtain various properties of the solution. The common problem is to find the least-norm solution: find $x$, such that<br /> $$\min_x \|x\|_2, \text{ such that } Ax = b.$$<br /> <br /> The solution is given by $x = A^\T (AA^\T)^{-1} b$, and it can be obtained by solving the $n \times n$ system $AA^\T y = b$ and computing $x = A^\T y$. <br /> The solution can also be written as $x = A^+ b$. <br /> <br /> To validate the correctness, consider any other $z$, such that $Az=b$. Then, $$<br /> (z-x)^\T x = (z-x)^\T A^\T (AA^\T)^{-1} b = (A(x-z))^\T (AA^\T)^{-1} b = (b-b)^\T (AA^\T)^{-1} b = 0,<br /> $$<br /> and consequently<br /> $$<br /> \|z\|^2 = \|z-x+x\|^2 = \|z-x\|^2 + 2(z-x)^\T x + \|x\|^2 = \|z-x\|^2 + \|x\|^2 \geq \|x\|^2,<br /> $$<br /> which shows that $x$ has smaller norm than any other solution $z$.<br /> <br /> Alternative derivation is to minimize $x^\T x$ with respect to $Ax=b$ using Lagrangian multipliers.<br /> <br /> = Weighted case =<br /> <br /> Weighted case is similar, except that the minimization problem becomes $$\min_x \|x\|_{2,w}, \text{ such that } Ax = b.$$<br /> The minimization can be written $\min_x \|Wx\|_{2}, \text{ such that } Ax = b$, or substituting $Wx = y$, we obtain the following<br /> ordinary least norm problem:<br /> $$\min_y \|y\|_{2}, \text{ such that } AW^{-1}y = b,$$<br /> with matrix $AW^{-1}$ and right hand side $b$, and with solution $x = W^{-1} y$.<br /> <br /> Thus, the final solution is given as $x = W^{-2} A^\T (A W^{-2} A^\T)^{-1} b$ and can be computed by first solving <br /> $n \times n$ system $AW^{-2}A^\T y = b$ and computing $x = W^{-2} A^\T y$.<br /> <br /> = Numerical solving =<br /> <br /> Techniques for numerical solving of (weighted) underdetermined systems are described below. As list of decompositions available in Eigen is available here:<br /> https://eigen.tuxfamily.org/dox/group__LeastSquares.html<br /> <br /> As when solving overdetermined systems if a lot of systems $Ax = b_i$ need to be solved, which only differ in the right hand side, the matrix $A$ can only be decomposed once, with<br /> subsequent solutions taking $O(n^2)$ time.<br /> <br /> == Normal equations ==<br /> <br /> We can solve the system directly, by solving the $n \times n$ system $AA^\T y = b$ using Cholesky, and computing $x = A^\T y$.<br /> <br /> This decomposition takes $\frac13 n^3 + 2n^2m$ time.<br /> <br /> == QR decomposition ==<br /> <br /> Writing $A^\T = Q_1R_1$ as when solving overdetermined systems, we can get the solution<br /> $$x = Q_1R_1 (R_1^\T Q_1^\T Q_1 R_1)^{-1}b = Q_1 R_1 R_1^{-1} R_1^{-\T} b = Q_1 R_1^{-\T} b.$$<br /> <br /> The computational costs are the same as in the overdetermined case.<br /> <br /> == SVD decomposition ==<br /> <br /> The pseudoinverse can be computed using SVD and the solution obtained as $x = A^{+} b$.<br /> <br /> The computational costs are the same as in the overdetermined case.<br /> <br /> == Decompositions ==<br /> <br /> Solving linear systems is not done by computing inverses of $A$ but by suitably decomposing $A$. <br /> Many decompositions are available, with $LU$ decomposition with partial pivoting being the go-to option.<br /> with $\frac23 n^3$ cost. Additionally, we have $QR$ and $SVD$ decompositions if more robust decompositions are desired, at a greater cost.<br /> <br /> See full the list of Eigen's available decompositions: https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html<br /> <br /> If many systems $Ax = b_i$ need to be saved with only different right hand sides, the decomposition of $A$ can be computed in $O(n^3)$ and<br /> stored with subsequent solutions only taking $O(n^2)$ time.<br /> <br /> If the matrix is symmetric and positive definite, Cholesky (called $VV^\T$ or $LL^\T$) decomposition can be used, with $\frac13 n^3$ cost.<br /> <br /> = Basics of solving overdetermined linear systems =<br /> == ordinary least squares ==<br /> <br /> Let $A$ be a matrix of size $n \times m$ and $b$ be a vector of size $n$, with $n \geq m$ and let $A$ have full rank. If $n=m$ the system $Ax=b$ is square, and can be treated as such, however the<br /> derivation below works as well.<br /> <br /> The system $Ax=b$ has $n$ equations for $m$ unknowns $x_1, \ldots, x_m$ and is overdetermined. The solution is most often sought in the least-squares sense, defining the solution $x$ to be the solution of the following minimisation problem<br /> $$ \min_x \|Ax-b\|_2. $$<br /> The solution can be obtained using calculus, by looking for stationary points of function $\|Ax-b\|_2^2 = (Ax-b)^\T(Ax-b)$ (which has the same minimum), and getting the system of<br /> normal equations $$A^\T A x = A^\T b.$$ The solution can be expressed theoretically as<br /> $x = (A^\T A)^{-1} A^\T b$.<br /> <br /> == Weighted least squares ==<br /> <br /> The weighted least squares problem has the same setup as the ordinary least squares, but we instead minimise the norm $\|Ax-b\|_{2,w}$, where $\|z\|_{2,w} = \sum_{i=1}^n (w_iz_i)^2$. We assume all weights $w_i$ are positive.<br /> <br /> By observing that $\|z\|_{2,w} = \|Wz\|_2$, where $W = \operatorname{diag}(w_1, \ldots, w_n)$, we can reduce the problem to ordinary least squares.<br /> Writing $\|Ax-b\|_{2,w} = \|W(Ax-b)\|_2 = \|WA x - Wb\|_2$, we obtain an ordinary least squares problem with matrix $WA$ and vector $Wb$, which has the <br /> normal equations $$A^\T W^2Ax = A^\T W^2b$$ and the solution $x = (A^\T W^2A)^{-1} A^\T W^2b$.<br /> <br /> = Numerical solving =<br /> <br /> Techniques for numerical solving of (weighted) least squares systems are described below. As list of decompositions available in Eigen is available here:<br /> https://eigen.tuxfamily.org/dox/group__LeastSquares.html<br /> <br /> As when solving linear systems if a lot of systems $Ax = b_i$ need to be solved, which only differ in the right hand side, the matrix $A$ can only be decomposed once, with<br /> subsequent solutions taking $O(n^2)$ time.<br /> <br /> == Normal equations ==<br /> The normal equations above can be directly used to solve the system. The matrix $A^\T A$ is symmetric and positive definite, and<br /> normal equations can be solved using Cholesky decomposition. However, this can be <br /> poorly behaved with respect to round-off errors since the condition number $\kappa(A^\T A)$ is the square<br /> of $\kappa(A)$.<br /> <br /> Complexity of Cholesky decomposition is $\frac{m^3}{3}$ and complexity of matrix multiplication $2m^2n$. To preform the Cholesky decomposition the rank of $A$ must be full.<br /> <br /> '''Pros:'''<br /> * simple to implement<br /> * low computational complexity<br /> <br /> '''Cons:'''<br /> * numerically unstable<br /> * full rank requirement<br /> <br /> == [https://en.wikipedia.org/wiki/QR_decomposition $QR$ Decomposition] ==<br /> The square $n\times m$ matrix $A$ is decomposed as <br /> $$<br /> A = QR = <br /> \begin{bmatrix} Q_1 &amp; Q_2 \end{bmatrix} <br /> \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1R_1,$$<br /> where <br /> $Q$ is a unitary matrix ($Q^{-1}=Q^T$) of size $n\times n$, $R$ trapezoidal and is the same size as $A$ ($n \times m$),<br /> $Q_1$ is also the same as $A$ ($n \times m$) and $R_1$ is a triangular $m \times m$ matrix.<br /> A useful property of a unitary matrices is that multiplying with them does not alter the (Euclidean) norm of a vector, i.e.,<br /> $$ \| Qx \| = \| x \|. $$<br /> Therefore we can say<br /> $$<br /> \| Ax- b \| = \| Q^\T (Ax-b) \| = \| Q^\T Ax - Q^\T b \|<br /> = \| Q^\T QRx - Q^\T b \| = \| [R_1; 0] x - [Q_1, Q_2]^\T b \|<br /> = \| R_1x - Q_1^\T b \| + \| Q_2^\T b \|<br /> $$<br /> Of the two terms, we have no control over the second, and we can render the first one<br /> zero by solving<br /> $$ R_1x = Q_1^\T b, $$<br /> which results in a minimum $x = R_1^{-1} Q_1^\T b$.<br /> <br /> When computing the $QR$ decomposition, only the ''thin'' decomposition $Q_1$, $R_1$ is needed to solve the least squares problem.<br /> <br /> The complexity of $QR$ decomposition is $2nm^2−2m^3/3$.<br /> <br /> &lt;strong&gt;Pros:&lt;/strong&gt; better stability in comparison to normal equations, &lt;strong&gt; cons: &lt;/strong&gt;higher complexity<br /> <br /> == [https://en.wikipedia.org/wiki/Singular_value_decomposition SVD decomposition] ==<br /> In linear algebra, the [https://en.wikipedia.org/wiki/Singular_value_decomposition singular value decomposition (SVD)]<br /> is a factorization of a real or complex matrix. <br /> <br /> The singular value decomposition of an $n \times m$ real or complex<br /> matrix $A$ is a factorization of the form $A= U\Sigma V^\T$, where<br /> $U$ is an $n \times n$ unitary matrix, $\Sigma$ is an $n \times m$<br /> rectangular diagonal matrix with non-negative real numbers on the diagonal, and<br /> $V$ is an $m \times m$ unitary matrix. The diagonal entries<br /> $\Sigma_{ii}$ are known as the singular values of $A$. The $n$ columns of<br /> $U$ and the $n$ columns of $V$ are called the left-singular vectors and<br /> right-singular vectors of $A$, respectively.<br /> <br /> The singular value decomposition and the eigendecomposition are closely<br /> related. Namely:<br /> <br /> * The left-singular vectors of $A$ are eigen vectors of $A A^\T$.<br /> * The right-singular vectors of $A$ are eigen vectors of $A^\T A$.<br /> * The non-zero singular values of $A$ (found on the diagonal entries of $\Sigma$) are the square roots of the non-zero eigenvalues of both $A^\T A$ and $A^\T A$.<br /> <br /> with SVD we can write $A$ as \[A= U\Sigma V^\T \] where are $U$ and $V$ again unitary matrices and $\Sigma$<br /> stands for diagonal matrix of singular values.<br /> <br /> SVD can be used to compute the pseudoinverse<br /> \[ A^{+} = \left( U\Sigma V^\T\right)^{-1} = V \Sigma^{-1}U^\T \]<br /> where $\Sigma^{-1}$ is trivial, just replace every non-zero diagonal entry by<br /> its reciprocal and transpose the resulting matrix. <br /> <br /> This can be used to solve the least squares problem, ordinary systems or underdetermined systems.<br /> Note that usually only the ''thin'' variant is needed.<br /> <br /> One can also choose a threshold $t$<br /> below which the singular value is considered as $0$ and truncate all <br /> singular values, whose magnitude is below $t$. Because we do not invert very small values,<br /> stability of the solution increases, at the cost of a slightly inexact solution.<br /> This approach is called &quot;truncated SVD decomposition&quot;, which is a form of regularization of badly conditioned problems.<br /> Another approach would be [https://en.wikipedia.org/wiki/Tikhonov_regularization Tikhonov regularization].<br /> <br /> SVD decomposition complexity \[ 2nm^2+2m^3 = O(n^3) \]<br /> <br /> &lt;strong&gt;Pros:&lt;/strong&gt; stable, &lt;strong&gt; cons: &lt;/strong&gt;high complexity</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=Ghost_nodes&diff=3554 2023-11-05T21:05:00Z

<p>Zigavaupotic: </p> <hr /> <div>Go back to [[Medusa#Examples|Examples]].<br /> <br /> See the [[Ghost nodes (theory)]] page for what ghost nodes and how they are used. We will use them in this example to reliably solve a 3D mixed Dirichlet and Neumann problem on an irregular domain.<br /> <br /> We will solve the problem<br /> <br /> $\nabla^2 u = 1 \text{ in } \Omega, \quad \frac{\partial u}{\partial n} = 0 \text{ on } \partial \Omega_{+}, \quad u = 0 \text{ on } \partial \Omega_{-}$<br /> <br /> where $\Omega = B(\boldsymbol{0}, 1) - B(\boldsymbol{1}, 1.5)$, $\partial \Omega_{+}$ is the part of the boundary with the nonnegative $x$ coordinate and $\partial \Omega_{-}$ the part of the boundary with negative $x$ coordinate.<br /> <br /> The code below constructs ghost nodes and remembers the corresponding ghost node for each boundary node.<br /> &lt;syntaxhighlight lang=&quot;cpp&quot;&gt;<br /> Range&lt;int&gt; gh = domain.addGhostNodes(dx, 0, neu);<br /> &lt;/syntaxhighlight&gt;<br /> It is equivalent to:<br /> &lt;syntaxhighlight lang=&quot;cpp&quot;&gt;<br /> Range&lt;int&gt; gh(domain.size(), -1); // Maps boundary nodes to their corresponding ghost nodes.<br /> for (int i : neu) { // Add ghost nodes for each Neumann boundary node<br /> gh[i] = domain.addNode(domain.pos(i) + dx*domain.normal(i), 0); // add node with type 0<br /> }<br /> &lt;/syntaxhighlight&gt;<br /> <br /> When assembling the system, we write two equations for each Neumann node. One is the Neumann BC and is written in the row corresponding to the boundary node. <br /> The other is the fact that the PDE also holds at the boundary node, but we write that equation in the row, corresponding to the ghost node. <br /> <br /> &lt;syntaxhighlight lang=&quot;cpp&quot;&gt;<br /> for (int i : interior) {<br /> op.lap(i) = 1.0; // Equation<br /> }<br /> for (int i : dir) {<br /> op.value(i) = 0.0; // Dirichlet<br /> }<br /> for (int i : neu) {<br /> op.neumann(i, domain.normal(i)) = 0.0; // Neumann<br /> // Equation holds also on the boundary, write it in the row corresponding to the ghost node.<br /> op.lap(i, gh[i]) = 1.0;<br /> }<br /> &lt;/syntaxhighlight&gt;<br /> <br /> The system is then solved as usual. The solution (without the ghost nodes) is shown below:<br /> <br /> [[File:solution_ghost_example.png|600px]]<br /> <br /> The code is available as [https://gitlab.com/e62Lab/medusa/blob/dev/examples/ghost_nodes/poisson_mixed_3D_ghost.cpp poisson_mixed_3D_ghost.cpp] in the Medusa examples.<br /> <br /> Go back to [[Medusa#Examples|Examples]].</div>

Zigavaupotic http://e6.ijs.si/medusa/wiki/index.php?title=How_to_build&diff=3553 2023-07-07T14:19:04Z

<p>Zigavaupotic: </p> <hr /> <div>=Installation=<br /> To make this work from plain Ubuntu installation, run<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> sudo apt-get install git g++ python3 cmake libhdf5-serial-dev doxygen graphviz<br /> git clone https://gitlab.com/e62Lab/medusa.git --branch dev --single-branch<br /> cd medusa<br /> ./run_tests.py -t<br /> &lt;/syntaxhighlight&gt;<br /> which installs dependencies, clones the repository, goes into the root folder of<br /> the repository and runs tests. This will build and run all tests. If this<br /> works, you are ready to go! Otherwise install any missing packages and if it<br /> still fails, raise an issue!<br /> Note: If you are told the packages cannot be located, try doing a sudo apt-get update.<br /> <br /> For instructions on how to use this library in you project, see<br /> [[Including this library in your project]].<br /> <br /> =Building=<br /> <br /> List of dependencies:<br /> <br /> * Build tools, like &lt;code&gt; cmake &gt;= 2.8.12&lt;/code&gt;, &lt;code&gt; g++ &gt;= 4.8&lt;/code&gt;, &lt;code&gt;make&lt;/code&gt;, &lt;code&gt;python3&lt;/code&gt;<br /> * [https://www.hdfgroup.org/ HDF5 library] for IO <br /> * &lt;code&gt;doxygen &gt;= 1.8.8 &lt;/code&gt; and Graphviz for generating the documentation<br /> <br /> Out of source builds are preferred. Run<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> mkdir -p build<br /> cd build<br /> cmake ..<br /> make<br /> &lt;/syntaxhighlight&gt;<br /> Note that you only have to run &lt;code&gt;cmake&lt;/code&gt; once, after that only &lt;code&gt;make&lt;/code&gt; is sufficient.<br /> <br /> Binaries are placed into &lt;code&gt;bin/&lt;/code&gt; folder. Tests can be run all at once via &lt;code&gt;make medusa_run_tests&lt;/code&gt; <br /> or individually via e. g. &lt;code&gt;make operators_run_tests&lt;/code&gt;.<br /> <br /> == Linker errors ==<br /> <br /> When trying out different classes, you might come along linker errors such as <br /> <br /> &lt;code&gt;<br /> Scanning dependencies of target cantilever_beam<br /> [100%] Building CXX object examples/linear_elasticity/CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o<br /> [100%] Linking CXX executable ../../../examples/linear_elasticity/cantilever_beam<br /> /usr/bin/ld: CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o: in function `main':<br /> cantilever_beam.cpp:(.text.startup+0x162): undefined reference to `void mm::FindBalancedSupport::operator()&lt;mm::DomainDiscretization&lt;Eigen::Matrix&lt;double, 2, 1, 0, 2, 1&gt; &gt; &gt;(mm::DomainDiscretization&lt;Eigen::Matrix&lt;double, 2, 1, 0, 2, 1&gt; &gt;&amp;) const'<br /> collect2: error: ld returned 1 exit status<br /> &lt;/code&gt;<br /> <br /> This is expected and is the result of some optimizations of compilation time. The medusa library can actually be included in two ways: as<br /> &lt;code&gt;#include &lt;medusa/Medusa_fwd.hpp&gt;&lt;/code&gt; or &lt;code&gt;#include &lt;medusa/Medusa.hpp&gt;&lt;/code&gt;. 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 &lt;i&gt;full&lt;/i&gt; Medusa library (the header without the &lt;code&gt;_fwd&lt;/code&gt;) 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 &lt;code&gt;medusa/bits/domains/FindBalancedSupport.hpp&lt;/code&gt;) and pay for whay you use. Or, add your instantiation among the already precompiled instantiations (located in &lt;code&gt;.cpp&lt;/code&gt; files, such as e.g. [https://gitlab.com/e62Lab/medusa/-/blob/dev/src/domains/DomainDiscretization.cpp this one]).<br /> <br /> == Building on macOS ==<br /> This method was tested on macOS Mojave 10.14.2.<br /> <br /> First install Xcode via App Store and then Xcode Command Line Tools with<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> xcode-select --install<br /> &lt;/syntaxhighlight&gt;<br /> <br /> After that, install all dependencies from homebrew<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> brew install cmake hdf5 boost python doxygen graphviz<br /> &lt;/syntaxhighlight&gt;<br /> <br /> Now you can clone and build the project using the following commands<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> git clone https://gitlab.com/e62Lab/medusa.git<br /> cd medusa<br /> mkdir build<br /> cd build<br /> cmake ..<br /> cd ..<br /> python3 run_tests.py -t<br /> &lt;/syntaxhighlight&gt;<br /> <br /> 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!<br /> <br /> ==HDF5==<br /> In order to use HDF5 IO you need the [https://www.hdfgroup.org/ HDF5 library].<br /> You can install it easily using the command &lt;code&gt;sudo apt-get install libhdf5-dev&lt;/code&gt;<br /> or &lt;code&gt;sudo pacman -S hdf5&lt;/code&gt;.<br /> <br /> Ubuntu places (at least on older versions) hdf5 headers and libraries in a weird folder <br /> &lt;code&gt;/usr/{lib, include}/x86_64-linux-gnu/hdf5/serial/&lt;/code&gt;. <br /> <br /> If you get an error like &lt;code&gt;HDF5 sample failed to compile. See errors above.&lt;/code&gt; during &lt;code&gt;cmake&lt;/code&gt; execution<br /> then the sample hdf test file located in &lt;code&gt;test/test_hdf_compile.cpp&lt;/code&gt; failed to compile. Perhaps it is good to make this file compile first, <br /> before tackling the whole project. <br /> <br /> If you get an error similar to &lt;code&gt;fatal error: hdf5.h: No such file or directory&lt;/code&gt;,<br /> then your compiler cannot see the HDF5 header files. Some distributions, notably (older) Ubuntu, place them into nonstandard folders<br /> &lt;code&gt;/usr/include/hdf5/serial/&lt;/code&gt; or &lt;code&gt;/usr/include/x86_64-linux-gnu/hdf5/serial/&lt;/code&gt;.<br /> Check these two folders or check your distributions hdf package for the locations of these files.<br /> After you have determined the location, add that directory to the include directories,<br /> using -I flag or in &lt;code&gt;CMakeLists.txt&lt;/code&gt; by using<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> include_directories(/usr/include/hdf5/serial/) # or your appropriate directory<br /> &lt;/syntaxhighlight&gt;<br /> If you wish to fix this problem permanently, you can create soft links to the headers in your &lt;code&gt;/usr/include&lt;/code&gt; directory, <br /> by typing<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> sudo ln -s /usr/include/hdf5/serial/* /usr/include<br /> &lt;/syntaxhighlight&gt;<br /> After this, there should be no compile time errors. If there are, please raise an issue.<br /> <br /> If you get error similar to &lt;code&gt;-lhdf5 not found&lt;/code&gt; and you have hdf5 installed, <br /> you might have to link the libraries into a discoverable place, like &lt;code&gt;/usr/lib/&lt;/code&gt; <br /> or add the above directory to the linker path. Similarly to above, check the &lt;code&gt;/usr/lib/x86_64-linux-gnu/hdf5/serial/&lt;/code&gt;<br /> directory and look for file &lt;code&gt;libhdf5.a&lt;/code&gt;. When found,<br /> specify the location using -L flag or &lt;code&gt;CMakeLists.txt&lt;/code&gt; by using<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> link_directories(/usr/lib/x86_64-linux-gnu/hdf5/serial/) # or your appropriate directory<br /> &lt;/syntaxhighlight&gt;<br /> or fix the problem permanently by soft-linking<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> sudo ln -s /usr/lib/x86_64-linux-gnu/hdf5/serial/* /usr/lib<br /> &lt;/syntaxhighlight&gt;<br /> <br /> == OpenMP ==<br /> Sometimes, OpenMP cmake errors might occure. This happens mainly due to limited multi-thread support. One can fix such issues, by adding the following code into their CMakeLists.txt<br /> <br /> &lt;syntaxhighlight lang=&quot;cmake&quot;&gt;<br /> find_package(OpenMP)<br /> if (OPENMP_FOUND)<br /> set (CMAKE_C_FLAGS &quot;${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}&quot;)<br /> set (CMAKE_CXX_FLAGS &quot;${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}&quot;)<br /> set (CMAKE_EXE_LINKER_FLAGS &quot;${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}&quot;)<br /> endif()<br /> &lt;/syntaxhighlight&gt;<br /> <br /> ==Linear Algebra==<br /> We use [http://eigen.tuxfamily.org/ Eigen] as our matrix library. See<br /> [http://eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html here] for use<br /> reference and documentation. For a quick transition from Matlab see<br /> [http://eigen.tuxfamily.org/dox/AsciiQuickReference.txt here].<br /> <br /> == Using Intel Math Kernel Library (MKL) ==<br /> Install [https://software.intel.com/en-us/mkl Intel MKL] and take not of installation directories.<br /> <br /> Proper include and link directories need to be set to use the Intel MKL. <br /> &lt;syntaxhighlight lang=&quot;cmake&quot;&gt;<br /> include_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/include) # change these to your installation path<br /> link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/lib/intel64)<br /> link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/lib/intel64)<br /> &lt;/syntaxhighlight&gt;<br /> <br /> Eigen has great support for MKL. You can see more detailed instructions [https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html on their website].<br /> To use MKL for math operations, define &lt;code&gt;EIGEN_USE_MKL_VML&lt;/code&gt; when compiling. You must also link<br /> the appropriate libraries and define &lt;code&gt;MKL_LP64&lt;/code&gt; for use on a 64bit system.<br /> With parallelism enabled, the configuration looks as follows<br /> <br /> &lt;syntaxhighlight lang=&quot;cmake&quot;&gt;<br /> target_compile_options(my_target PRIVATE &quot;-fopenmp&quot;)<br /> target_compile_definitions(my_target PUBLIC EIGEN_USE_MKL_VML MKL_LP64)<br /> target_link_libraries(my_target medusa mkl_intel_lp64 mkl_intel_thread mkl_core pthread iomp5)<br /> &lt;/syntaxhighlight&gt;<br /> <br /> 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].<br /> <br /> == Using Intel C/C++ Compiler ==<br /> <br /> In order to use Intel's compiler when compiling Medusa, you have several standard optionas for &lt;code&gt;cmake&lt;/code&gt;.<br /> Make sure compilers and installed and in your &lt;code&gt;PATH&lt;/code&gt; by running &lt;code&gt;which icc&lt;/code&gt;, which <br /> should return something like &lt;code&gt;/opt/intel/bin/icc&lt;/code&gt;.<br /> <br /> You can define the compiler when *first* calling cmake like so<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> cmake .. -DCMAKE_C_COMPILER=$(which icc) -DCMAKE_CXX_COMPILER=$(which icpc) <br /> &lt;/syntaxhighlight&gt;<br /> If this is not your first call, remove the &lt;code&gt;build&lt;/code&gt; directory and start anew.<br /> <br /> You can also set the &lt;code&gt;CXX&lt;/code&gt; and &lt;code&gt;CC&lt;/code&gt; bash variables. Before calling<br /> &lt;code&gt; cmake &lt;/code&gt; for the first time you have to export the following:<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> export CXX=&quot;icpc&quot;<br /> export CC=&quot;icc&quot;<br /> &lt;/syntaxhighlight&gt;<br /> <br /> &lt;!--<br /> You can also compile it directly for Intel® Xeon Phi™ Coprocessor. You do this by adding &lt;code&gt;-Dmmic=ON&lt;/code&gt;<br /> flag to the &lt;code&gt;cmake&lt;/code&gt; command:<br /> &lt;syntaxhighlight lang=&quot;bash&quot;&gt;<br /> cmake .. -Dmmic=ON -DCMAKE_C_COMPILER=$(which icc) -DCMAKE_CXX_COMPILER=$(which icpc) <br /> &lt;/syntaxhighlight&gt;<br /> <br /> <br /> &lt;b&gt;Note:&lt;/b&gt; All features that depend on system third-party libraries are not available on MIC (Many Integrated Core).<br /> This includes:<br /> <br /> * HDF class in &lt;code&gt;io.hpp&lt;/code&gt;<br /> <br /> --&gt;</div>

Zigavaupotic