Editing Appendix/Ramblings/ForCohlHoward (section)

=Considering Utility of Finite-Element Techniques=

==Joel's Opening Comments==
Among the introductory chapters of our [[H_BookTiledMenu|MediaWiki-based discussion of ''Self-Gravitating Fluids'']], you will find a formal [[PGE|list of the relevant "Principal Governing Equations]]" and an accompanying discussion of, for example, the [[PGE/Euler|Euler equation]].  Here, we pull from that discussion the so-called,

<div align="center">
<span id="ConservingMomentum:Conservative"><font color="#770000">'''Conservative Form'''</font></span><br />
of the Euler Equation,

<math>\frac{\partial(\rho\vec{v})}{\partial t} + \underbrace{\nabla\cdot [(\rho\vec{v})\vec{v}]}_{\mathrm{advection}} 
= 
\overbrace{- \nabla P - \rho \nabla \Phi}^{\mathrm{source}}</math>

[<b>[[Appendix/References#BLRY07|<font color="red">BLRY07</font>]]</b>], p. 8, Eq. (1.31)
</div>

which contains the divergence of the "[https://en.wikipedia.org/wiki/Dyadics dyadic product]" or "[https://en.wikipedia.org/wiki/Outer_product outer product]" of the vector momentum density and the velocity vector, and is sometimes written as, <math>\nabla\cdot [(\rho \vec{v}) \otimes \vec{v}]</math>.  How do we rewrite this ''continuum'' equation in a ''discrete'' form that will allow us (using digital computers) to accurately integrate this equation &#8212; in tandem with the other "Principal Governing Equations" &#8212; to determine how the fluid momentum at every point in space evolves in time?  Answering this question is a particularly challenging task because the so-called ''advection'' term on the left-hand-side is inherently nonlinear in the velocity.

Over the years, applied mathematicians &#8212; including, for example, engineers, physicists, astronomers, aerodynamicists &#8212; have devised a large array of schemes, the most popular falling under one of the following three broad (Wikipedia-acknowledged!) categories:
<ul>
  <li>
[https://en.wikipedia.org/wiki/Finite_difference_method Finite-Difference Method] (FDM)
  </li>
  <li>
[https://en.wikipedia.org/wiki/Finite_volume_method Finite-Volume Method] (FVM) &nbsp; &nbsp; <math>\Leftarrow</math> &nbsp; Our category of choice, historically (see immediately below)
  </li>
  <li>
[https://en.wikipedia.org/wiki/Finite_element_method Finite-Element Method] (FEM)
  </li>
</ul>

When I googled "finite element vs finite difference " on 13 October 2022, I was pointed to, for example,
<ul>
  <li>
[https://www.sciencedirect.com/science/article/pii/B9780122532504500191 Comparison of Finite-Element and Finite-Difference Methods]
  </li>
  <li>
[https://www.machinedesign.com/3d-printing-cad/fea-and-simulation/article/21832072/whats-the-difference-between-fem-fdm-and-fvm What's the Difference Between FEM, FDM, and FVM?]
  </li>
</ul>

===Our Choice, Historically===

Over the past four decades, the nonrelativistic fluid simulations that have been performed by LSU's astrophysics group have been carried out using various renditions of an ''explicit'' <b>Finite-Volume Method</b>.  The explicit FVM method was chosen primarily because &hellip;
<ol>
  <li>
It was the method to which I was introduced by my [[Appendix/Ramblings/MyDoctoralStudents#Years_1976_-_1978|doctoral dissertation advisors]], Peter Bodenheimer and David Black, that was considered appropriate to the type of fluid-flows problems that were the focus of my dissertation research.  Note that Bodenheimer is the first author of the book that is [[#Joel%27s_Opening_Comments|referenced above]] as [[Appendix/References#BLRY07|<font color="red">BLRY07</font>]].
  </li>
  <li>
Time-integration is fairly straightforward because all terms (except the time derivative) are specified entirely in terms of the current (as opposed to advanced) time.
  </li>
  <li>
The nonlinear ''advection'' term on the LHS is written in a way that guarantees global momentum conservation, if the ''source'' terms on the RHS of the Euler equation sum to zero.
  </li>
  <li>
The amount of computer time that is required to take each step forward in time is typically much smaller when using an ''explicit'' as opposed to ''implicit'' integration scheme.
  </li>
</ol>

A major disadvantage of our decades-old, explicit-FVM approach is that, in order for the scheme to be numerically stable, roughly speaking each integration time-step <math>\Delta t</math> must be smaller than the smallest value of the ratio, <math>\Delta x/c_s</math>, as determined across the entire grid, where <math>\Delta x</math> is the grid spacing and <math>c_s</math> is the fluid sound-speed.  (Note that <math>\Delta t</math> is essentially the sound-crossing time across the smallest grid cell.)  This constraint on the integration time-step is computationally burdensome when attempting to employ a high-resolution (small <math>\Delta x</math>) grid and/or when attempting to model environments where the fluid is relatively hot and therefore characterized by a high sound speed.  Even worse, <b>it is ''impossible'' to model the evolution of ''incompressible'' fluids using an explicit-FVM scheme</b> because, in effect in such fluids, sound travels at an infinite speed; that is, <math>c_s = \infty</math> so <math>\Delta t</math> goes to zero.

===Quandary===

Because we have chosen over the past four decades to employ an explicit-FVM approach in our examination of the stability of rotating, self-gravitating fluids, we have been constrained to examine the behavior of ''compressible'' rather than ''incompressible'' fluids.  This is perfectly acceptable &#8212; and even desirable &#8212; in the context of studies of star-formation, for example, because we are confident that star-forming gas clouds and young protostars are composed of fluids that obey compressible equations of state.  But it is difficult to critique the validity &#8212; and therefore relevance &#8212; of the results of such compressible-fluid simulations because <b>the published literature contains virtually no linear-stability analyses of compressible fluid systems</b> against which our nonlinear simulations can be compared.

On the other hand, as Chandrasekhar details in [[Appendix/References#EFE|<font color="red">EFE</font>]], the literature is packed full of studies that have quantitatively assessed the relative stability of rotating, self-gravitating, ''incompressible'' fluid systems.  (And, associated with the arrival of Stephen Sorokanich at NIST, we expect this literature to expand significantly.)  So, the validity &#8212; and relevance &#8212; of our numerical simulations would be on a much firmer foundation if we had the ability to model the evolution of incompressible fluid systems.  As an additional reward, a numerical simulation tool of this type would allow us to follow the ''nonlinear'' development of instabilities (in all Riemann ellipsoids, for example) whose onset has been predicted only from ''linear'' stability analyses.

===Perhaps We Should Consider the FEM===

Although I have had virtually no experience developing or using numerical algorithms that fall into the category of the Finite-Element Method (FEM), it is my understanding that such algorithms have the following features:
<ul>
  <li>
They can be used to follow the evolution of ''incompressible'' &#8212; and, for example, uniform-density &#8212; fluid systems.
  </li>
  <li>
They generally (always?) incorporate an ''implicit'' time-integration scheme in which case time steps are not constrained by &#8212; and can be much larger than &#8212; the often severely limiting sound-crossing times encountered in ''explicit'' FVM schemes.  
  </li>
  <li>
They are particularly good at identifying ''surfaces'' and modeling the time-evolution of surface distortions.
  </li>
</ul>

If you (Howard) and Sorokanich decide to build a hydro-code based on the Finite-Element Method, <font color="red">a potentially helpful reference</font> is:  
<div align="center">
{{ Meier99figure }}
</div>
I interacted with David Meier at a couple of different meetings back at the time he was developing this FEM algorithm. I'm not sure what astrophysics problems he tackled with this code or, ultimately, what he did with it; for example, I have never seen "the second paper in this series" to which he refers in the abstract.

In the context of our examination of the stability of Riemann ellipsoids, an FEM might permit us to reduce the dimensionality of the problem.  For example, ignore details associated with the interior of the 3D configuration and focus on modeling the distortion of its 2D surface.  

We might even be able to model (for the first time!) the topological transformation that is experienced by the surface when a rapidly spinning, incompressible ellipsoid ''fissions'' into a pair of disconnected tear-drop shaped objects in orbit about one another.  See, for example, the dissertation research of [https://www.semanticscholar.org/paper/A-hybrid-variational-level-set-approach-to-handle-Walker/7ee624ab9ffe45241cb0e5a0ce0898a6da201e7b Shawn W. Walker (2007)] (LSU Mathematics &amp; CCT), or more recent papers that Walker has coauthored with his dissertation advisor, Ricardo H. Nochetto (U. Maryland, Mathematics).

==Question Regarding Pressure-Poisson Equation==

<table border="1" align="center" width="80%" cellpadding="10"><tr><td align="left">
In an email dated 07 January 2023, Howard asked, "Do you think the pressure Poisson formulation of the Euler equations coupled with a 3D Poisson solver for gravity and evolution of the momentum equation might be our best strategy to evolve unstable Riemann S-type ellipsoids? We’ve already identified two different pieces of software (Fenics and COMSOL) which look like they might have this ability using finite-element. We also now have access to software which should be able to build good grids for FEM. Thanks! "
</td></tr></table>

Joel's response (11 January 2023) &hellip;

Your suggested "coupled Poisson formulation" sounds like a smart way to proceed.  Note, however, that I say this with limited authority; I have never written a code to solve the Euler equations that is based on a pressure-Poisson formulation.  My confidence in offering support for this approach comes from the following historical note.


In the early 1990s, in collaboration with Priya Vashishta, I helped the LSU Physics & Astronomy department secure approximately $0.8 million from the Louisiana Board of Regents to purchase [[Appendix/Ramblings/MyDoctoralStudents#Years_1988_-_1994|an 8K-node SIMD architecture MasPar MP1]].  As you know, our astrophysics group developed a 3D Poisson solver that was extraordinally well suited to the 8K-node hardware layout of the MasPar.  Around the same time that we received this ''state'' funding, I heard that the U.S. Department of Energy (DOE) was expecting to fund a limited number of large-scale fluid-flow (and heat-transfer) simulations in support of the aircraft industry's efforts to design better turbine engines.  LSU's research office organized a meeting of (primarily engineering) faculty to see if a group could be put together to submit a viable research proposal to the DOE.  I attended this meeting, just for grins; but it turned out to be a very smart decision.  At the meeting, I met [https://www.lsu.edu/eng/mie/people/formerfaculty/acharya.php Sumanta Acharya] (and [https://www.lsu.edu/eng/mie/people/faculty/nikitopoulos.php Dimitris Nikitopoulos]) and learned about his multidimensional simulations of incompressible fluids; he was looking for funding opportunities that would support his efforts to scale from 2D to 3D simulations.  I also learned that he usually solved the Euler equations using a pressure-Poisson approach.  We both realized that my group's successful implementation of a 3D Poisson solver (for gravity) on the MasPar could serve as an illustration of how his group could likely efficiently solve the pressure-Poisson equation in 3D.  We wrote the proposal to DOE, and it was funded at $0.5 million.  Thus began a very successful, multi-decade collaboration between my group and Sumanta's mechanical engineering group.


So &hellip; I am fairly confident that the approach you take to solve the gravitational-Poisson equation will help you understand how to solve the pressure-Poisson equation (or vise versa) on the finite-element grid.  Note, however, that to my knowledge, Sumanta's simulations always involved a finite-difference, rather than a finite-element, methodology.