Stability Analysis Experiments

[lukem12345]

This issue is for tracking the results of running pressure simulations on spherical meshes.
For each trial run, we should record the parameters and mesh used, as well as whether the solver was stable, and if any artifacts were observed in the simulation results (if the simulation was run).

[JarodGuzik]

mu=0, alpha=0.5, beta=20, gamma=20

[lukem12345]

From discussion with @bosonbaas ,

We should:
Add a drag term, and look at its steady state of velocity over the initial pressure field. This can be done experimentally, no changes to Decapodes.jl are necessary for this.
(A more physically realistic scenario would be to have the pressure field change over time, but this would require some larger changes.)

Negative pressure gradients do not result in unstable physics. Is this okay?
Yes.
Having a negative acceleration just changes your interpretation of high vs. low pressure.

There is a choice to have gradient either point from high to low pressure. The convention is to have low-to-high, and this is what we have in the code now. Confusion on this would result in scalars with the wrong sign.

How do we handle plotting of 1Forms?
Typically, we would take the 1forms, and smooth it into a continuous vector field and plot that.
The sharp operator might work out of the box. It would give us a vector on (likely dual?) points.
Hirani offers about 8 possible definitions of sharp. We believe the one we chose should do what we want here.

In matlab, we often have 3 separate plots, one for x component, one for y, and one for z.
We want plots in this case to be for theta\hat and phi\hat of the sharpened 1Forms.

When running experiments, we have to change alpha every time we change the mesh, since it is affected by lengths of edges.
Does the d_0 operator automatically accommodate changes in mesh size?
Yes, all similar operators are completely mesh independent.

Observe that since Velocity is a primal 1Form, the affect of pressure on velocity should also be a 1Form.
So if pressure is a primal 0Form, you would just take d of it, (and scale).

We note that normally in physics, the Laplacian applied to velocity gives you something in 1/(ms) units. In the DEC we mistakenly considered the Laplacian applied to velocity to be something that you can add to velocity (i.e. is also in m/s), but our units of pressure and velocity were not consistent.
Velocity in the DEC should actually be modeled in units of meters squared per second.

Using Unitful was helpful for debugging purposes. But they ran into some trouble with implicitly defined meshes.
Whenever you provide type info in your EmbeddedDeltaSet2D, you do e.g. m{Float64} (or whichever the appropriate syntax is) in place of just Float64.

A test that the scalars scale correctly for each icosphere subdivision would be to run the same physics on two different icospheres and confirm the simulation is the same (up to coarseness).

[JDCampII]

Ran three simulations using different mesh subdivisions. All had mu = 0, alpha = 0.5, beta = 100, gamma = 100.

Results:

• 3 Subdivisions, Simulation Time of ~45 secs, No Error/Warning Messages

• 4 Subdivisions, Simulation Time of ~12 mins, No Error/Warning Messages

• 5 Subdivisions, Simulation Time of ~1 hr, 39mins, Following Warning Message: Warning: dt(3.552713678800501e-15) <= dtmin(3.552713678800501e-15) at t=1.2868500015237205. Aborting. There is either an error in your model specification or the true solution is unstable.

[JDCampII]

Overall Stability Results:

• $0<\alpha<\infty(?)$
• $0<\beta<1e7$
• $0<\gamma<2$

Stability Results:

• $\alpha = \gamma = 0, \beta = 10e8$ - Simulation took too long to reasonably solve. Assumed unstable.
• $\alpha = \gamma = 0, \beta = 100$ - Stable. Outcome was a static spot on sphere.

• $\alpha = \gamma = 0, \beta = 10000$ - Stable. Extremely slow diffusion from initial source.
• $\alpha = \gamma = 0, \beta = 1e6$ - Stable. Fast diffusion with "blip" afterwards on the left-side of sphere.

• $\alpha = \gamma = 0, \beta = 1e8$ - Took too long to solve. Assumed unstable.
• $\alpha = \gamma = 0, \beta = 1e7$ - Stable. Near instantaneous diffusion. Long solving time.
• $\alpha = \gamma = 0, \beta = 5e7$ - Took too long to solve. Assumed unstable.
• $\alpha = \beta = 0, \gamma = 100$ - Unstable. Model self-deletes (sphere destroys itself).
• $\alpha = \beta = 0, \gamma = 10$ - Unstable. Mesh tends to $-\infty/\infty$

• $\alpha = \beta = 0, \gamma = 1$ - Stable. Diffuses normally.

• $\alpha = \beta = 0, \gamma = 3$ - Stable until diffusion reaches the poles. Then the model begins to tend towards $-\infty/\infty$.

[JarodGuzik]

Run time to solve with 4 subdivisions was 3.5 minutes
constants = (k=1, μ=-.001, α=0.5, βc=1e4, γc=1.0, βₚ=1e4, γₚ=1.0)

[JarodGuzik]

5 subdivisions still appears to be unstable. When I run with the same constants as above, I get this error
Warning: No strict ticks found
└ @ PlotUtils C:\Users\guzikjar.julia\packages\PlotUtils\igbcf\src\ticks.jl:191
ERROR: InexactError: trunc(Int64, NaN)

[Snider-wc]

my_consts = (k=1,μ=-10,α=1,βc=1e4,γc=1,βₚ=1e4,γₚ=1)
meshgrid division set at: 4

[lukem12345]

Filename: weather_alpha_1000_11_29_b_1e5_gc_m1_gp_m1_equator

(k=1,μ=-1,α=1000,βc=1e5,γc=-1,βₚ=1e5,γₚ=-1)

Setting the gammas to -10 will speed up the physics.

[lukem12345]

Commit 758296 performs a multi-species simulation, both of which governed individually by Navier Stokes. These species do not interact, so this essentially running two single-species simulations at once.

Interaction of species will properly be handled by the collision terms from Schunk, which mostly consist of linear sums of differences of velocity 1-forms. Although interaction of species can of course be encoded by a multitude of less physically-rigorous strategies for development purposes.

[jpfairbanks]

This looks really good. The different species should have different material properties (N2 and O2) for example so we should be able to run this simulation of two independent species and see that pressure waves propagate slower in different media based on the speed of sound. Then we can add the collision terms. I think that the sum over differences of one forms approach based on shunck seems reasonable to me. It would be nice to have a slick DEC way of framing that, but I think that just a sum over differences fine for now.

Once we have multiple species fluids going, we can add the E&M terms for charged particles.