Week of 2025.3.10¶

Last week I tried to get the slip PR in a reviewable and mergable state. One of the points that Eric made was that I should make sure to update the get_jacobian function with my new state variable.

Related to the above, the full jacobian should be a bit larger when adding another dependent variable (slip). I doubt the jacobian function would modified by Fred to > account for this. In fact, I'm not sure if he broke the ability to do operator splitting and implicit time stepping when he added slip as a new variable. The truly correct thing would be to make sure the jacobian is correct, but that could be a lot of work, especially in the super general framework where Fred allows for nonlocal dependence on fields in the term that is treated implicitly. I'm not sure that it's worth your time to figure this out and implement the needed changes. But at least we ought to include some kind of comment in the code and/or repo about this. Personally, I would recommend that you not update this full jacobian, but instead make note of it. Then as part of the new set of code changes, add support for an implicit term that is local in space. That set of changes will include the correct treatment of slip-dependence of the implicit term.

After a bit of investigation at the end of last week, I found that in order to handle one of the source terms implicity it was nessesary to mark a specific source as implicit and also change the timestepping method to stang, which is basically a way to solve subproblems of a differential equation seperately and still guarantee second order accuracy.

Before we merge the slip code, we want to make sure the general jacobian is correct.

However, the more important goal is to implement "local" implicit updates rather.

1.0 Imports¶

In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt

# Modify base path for depending on your file structure.
BASE_PATH = "/Users/paxton/git"

# Specify path where .pkl files are located
target_dir = f"{BASE_PATH}/quail_volcano/scenarios/simple_1D_test"
# Specify path for Quail source code
source_dir = f"{BASE_PATH}/quail_volcano/src"
# Change to working directory
os.chdir(target_dir)


# Import quail modules
os.chdir(source_dir)

import meshing.tools as mesh_tools

import numerics.helpers.helpers as helpers
import numerics.timestepping.tools as stepper_tools

import physics.zerodimensional.zerodimensional as zerod
import physics.euler.euler as euler
import physics.navierstokes.navierstokes as navierstokes
import physics.scalar.scalar as scalar
import physics.chemistry.chemistry as chemistry
import physics.multiphasevpT.multiphasevpT as multiphasevpT

import processing.readwritedatafiles as readwritedatafiles
import processing.post as post
import processing.plot as plot
import processing.animate as animate

import solver.DG as DG
import solver.ADERDG as ADERDG
import solver.tools as solver_tools

import time
from IPython.display import HTML
import multiprocessing as mp  
from multidomain import Domain, Observer

os.chdir(target_dir)

1.1 Getting the slip code to a place where it is merge ready¶

  • When the default viscosity is set, set all partial derivatives with respect to viscosity equal to zero.
  • I did a couple code reviews and merged.
  • I am still a little worried I messed up the jacobian for the implicit assumption.
  • In the future, I should also add better regression testing.

2.1 Adding a jacobian for conduit wall drag¶

In the future, in order to support implicit timestepping for the conduit wall friction source term I will want to add a get_jacobian for that source. The derivate of momentum with respect to slip aught to look something like:

$$ \begin{align} \frac{d_{mom}}{d slip} =& \frac{d \tau}{d slip} \frac{d mom}{d \tau} \\ \tau(slip) =& \tau_r - (\tau_r - \tau_p) e^{-slip/D_c} \\ momentum =& \frac{d \rho u}{d t} = - \frac{2 \tau}{R} \\ \frac{d \tau}{d slip} =& - \frac{1}{D_c} (- (\tau_r - \tau_p) e^{-slip/D_c}) \\ \frac{d_{mom}}{d \tau} =& - \frac{2}{R} \\ \frac{d_{mom}}{d slip} =& - \frac{1}{D_c} \frac{2}{R} (\tau_r - \tau_p) e^{-slip/D_c} \end{align} $$

When we are using the linear model for slip, we the derivative becomes a piecewise function: $\frac{d \tau}{d_{slip}} = \frac{\tau_r - \tau_p}{D_c}$ for slip < $D_c$ and 0 for all other values.

3.1 Verify that I did not break anything really obvious with my changes¶

In [3]:
ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_dynamic")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test
Out[3]:
Your browser does not support the video tag.

4.1 Review of timesteps and integrators¶

I found some nice review videos from Steve Brunton which review FE, BE, error analysis, and Runge-Kutta integrators.

Some relevant terms:

  • stiff: unstable except for very small time steps.
  • explicit numerical method: Calculates $X_{K+1}$ as a function of $X_k$, for instance forward euler: $X_{k+1} = X_k + \Delta t f(X_k)$.
  • implicit numerical method: Calculates $K_{k+1}$ as a function of $X_{k+1}$, for instance backward euler: $X_{K+1} = X_k + \Delta t f(X_{k+1})$

Generally, implicit methods are computationally more complex but have much better stability conditions.

5.1 Cleanup tasks¶

  1. Add flag to switch between exponential and linear model for tau.
  2. Help Mario get up and running with slip for his animation.

6.1 Review parameters now that the code is all merged.¶

$$ \begin{align} M \ddot{s} = A \Delta p(s(t)) + 2 \pi R L (\tau_p - \tau(s)) \\ M \ddot{s} = \frac{-A^2s(t)}{S} + 2 \pi R L (\tau_p - \tau(s)) \end{align} $$

Let's consider the two primary competing forces for stability. First there is the shear stress drag as a result of slip. That can be express as:

$$ \begin{align} F_d &= 2 \pi R L (\tau_p - \tau(s)) \\ \end{align} $$

Note that $\tau(s) = \tau_r - (\tau_r - \tau_p) \exp{-s(t)/D_c}$. Also note that with a first order taylor expansion, $\exp(-s(t)/D_c) = 1 - s(t)/D_c$. So a first order estimate of $\tau(s)$ is: $\tau(s) = \tau_p + (\tau_r - \tau_p)(\frac{s(t)}{D_c})$.

Applying that approximation we get:

$$ \begin{align} F_d &= 2 \pi R L (\tau_p - \tau_r)(\frac{s(t)}{D_c}) \\ \end{align} $$

Now, let's consider the force applied by the pressure on the left hand side of the domain

$F_p = \frac{-A^2 s(t)}{S}$ where $S = \frac{A L_{magma}}{K}$ so that becomes:

$$F_p = \frac{-A K s(t)}{L_{magma}}$$

Now, to think about stability lets call the system unstable when:

$$ \frac{\frac{d F_d}{d s}}{\frac{d F_p}{d s}} < 1$$

The derivatives are easy to take:

$$ \begin{align} \frac{d F_d}{d s} &= 2 \pi R L_{plug} (\frac{\tau_p - \tau_r}{D_c}) \\ \frac{d F_p}{d s} &= \frac{A K}{L_{magma}} \end{align} $$

$$ \begin{align} \frac{2 \pi R L_{plug} (\frac{\tau_p - \tau_r}{D_c})}{\frac{A K}{L_{magma}}} > 1 \\ \frac{2 \pi R L_{plug} L_{magma} (\tau_p - \tau_r)}{A K D_c} > 1 \\ \end{align} $$

So we conclude that the system is become more unstable when slip is the range $(0, D_c)$ and

$$ D_c < \frac{2 L_{magma} L_{plug} (\tau_p - \tau_r)}{R K} $$

Plugging in our numbers we get:

$$ D_c < \frac{2 * 100m[m] * 900 [m] * (3e5 [Pa] - 1e5 [Pa])}{10 [m] * 1e9 [Pa]} $$

$$ D_c < 3.6m $$

What does this tell us? It says something like: once we start slipping, if $D_c < 3.6$ we will continue to slip at least 3.6m until slip > D_c. At that point we will be at constant tau and the system should stop slipping.

In the below example, the system only slips to 2m despite our analytical analysis predicting 3.6m. I think this may be an issue with our frictional force not capping at the existing momentum which causes the system to go negative.

In [6]:
ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_dynamic")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test
Out[6]:
Your browser does not support the video tag.

7.1 How to prevent frictional force from being excessive¶

7.1.1. Naive solution is just to cap the friction at the amount of momentum avalalible.¶

S[:, :, physics.get_momentum_slice()] = - np.minimum(volumetric_friction, S[:, :, physics.get_momentum_slice()])

Okay... I tried this before it results in a really screwy solution.

In [8]:
ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_v1")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test
Out[8]:
Your browser does not support the video tag.

7.1.2 Use an implicit scheme + Rewrite tau so that it is dependendant on velocity.¶

Basically, write a step function.