Week of 2025.2.3 Experiments

Fred added a new advection state variable to the quail repo in two commits.

1. The first commit implements the new conserved variable, and set up a simulation with the new variable in a 1D domain with rigid walls on both ends, and no source terms.
2. The second commit implements a new placeholder source term that does nothing (except calculating pressure, temperature, etc. and making them available for further calculations), and source terms to the input file

As a first step, I worked off Fred's branch to set up a steady state flow simulation where I was able to calculate the velocity of the flow based on conservation of momentum.

As a next step, I am going to begin manipulating the new advection quantity and use the method of characteristics to verify that the new quantity is correctly implemented.

First, let's make sure all our imports are set up so we can plot the steady state flow example I was previously working with.

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. Review the current state of our system

To start out with, I update the animation of our system to include the new state variable arhoX implemented by Fred. We can see it is encoded as pDensityX in our state indices.

solver = readwritedatafiles.read_data_file(f"steady_state_flow/test_output_0.pkl")
physics = solver.physics
physics.state_indices

{'pDensityA': 0,
 'pDensityWv': 1,
 'pDensityM': 2,
 'XMomentum': 3,
 'Energy': 4,
 'pDensityWt': 5,
 'pDensityC': 6,
 'pDensityFm': 7,
 'pDensityX': 8}

In the original configuration, the arhoX term is $0$ for all $(x, t)$.

Let's also start plotting every relevant quantity: pressure, visocosity, velocity along with the all the states above.

These graphs below define IC of the new state as $0$, leave boundary conditions undefined, and define no source term.

ani = animate.animate_conduit_pressure("steady_state_flow", iterations=100, file_prefix="test_output")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

2. Sinusodial initial conditions

Next, lets run that the same simulation, but with sinusodial initial conditions. Specifically:

Uq[:, :, iarhoX] = arhoX * (1+ - np.sin(2* np.pi* x/1000)

ani = animate.animate_conduit_pressure("steady_state_flow", iterations=100, file_prefix="sinusoidal_IC_test")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

2.1 Characteristic curves

I found this pretty good video to review the method of characteristics. The advection equation can be written as:

$$ \frac{\partial q}{\partial t} + u \frac{\partial q}{\partial x} = S$$

Applying the method of characteristics, we observe that:

$\frac{\partial q}{\partial t} = S$ along the curve $\frac{dx}{dt} = u$

That implies that our solution is of the form:

$$ q(x, t) = st + f(x-ut) $$

def plot_characteristics(folder, file_prefix, u=2, N_m = 200, N_t = 100, T_end = 40, X_start=-1000):

    X = np.zeros((N_t, N_m)) 
    F = np.zeros((N_t, N_m)) 

    for i in range(0, N_t):
        solver = readwritedatafiles.read_data_file(f"{folder}/{file_prefix}_{i}.pkl")
        arhoX = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityX")]
        f, _, = physics.get_conv_flux_interior(solver.state_coeffs)

        X[i, :] = arhoX.ravel()
        F[i, :] = f[:,0,8,0].ravel()


    fig, ax = plt.subplots(figsize=(6,6))
    im = ax.imshow(X, extent=[-1000, 0, T_end, 0], aspect='auto')

    fig.colorbar(im)

    u = 2

    # Plot the charecteristics we expect the solution to be constant along. 
    plt.axline((-1000, 0), slope=1/u, color="black", linestyle=(0, (5, 5)))
    plt.axline((-800, 0), slope=1/u, color="black", linestyle=(0, (5, 5)))
    plt.axline((-600, 0), slope=1/u, color="black", linestyle=(0, (5, 5)))
    plt.axline((-400, 0), slope=1/u, color="black", linestyle=(0, (5, 5)))
    plt.axline((-200, 0), slope=1/u, color="black", linestyle=(0, (5, 5)))

    ax.set_ylabel('Time (s)')
    ax.set_xlabel('Position (m)')
    ax.invert_yaxis()
    ax.set_title("New state variable (pDensityX) over time and space")

folder = "steady_state_flow"
file_prefix = "sinusoidal_IC_test"

plot_characteristics(folder, file_prefix)

We can see that the there is a lot of diffusion happening which makes it a little hard to tell if the initial conditions remain constant along the characteristic curves. Here is what Fred said about the diffusion:

The diffusion/smoothing you see from the sinusoidal profile is from numerical error. Since we're using SolutionOrder of 0, the numerical error is primarily diffusive, and smooths out rough edges over time. This numerical error decreases as the grid is refined and the timestep is made smaller, but is something we have to live with. By going to higher SolutionOrder the diffusive error can be improved as well. Higher solution order corresponds to using higher-order polynomials to represent the solution within each cell.

I so far have been unable to increase the order of the solution, but I am trying a simulation with more spatial and temporal steps to see if that improves the situation at all. Let's return to this conversation about numerical errors later. For now, lets continue to experiment with the source term and boundary conditions.

3. Time varying boundary conditions.

For the moment, I am going to set the boundary condition only on one side:

UqB[:,:,8] = 2*np.sin(2 * np. pi *t/10)

I set the IC equal to arhoX = 1 at all positions, and I continued with no source term.

ani = animate.animate_conduit_pressure("steady_state_flow", iterations=100, file_prefix="time_varying_bc_test")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

folder = "steady_state_flow"
file_prefix = "time_varying_bc_test"

plot_characteristics(folder, file_prefix)

4. Source terms

Now, that we have experimented with both the ICs and the BCs, lets experiment with a source terms. First, we will have a constant source term of arhoX = 0.05

4.1 Constant source term

ani = animate.animate_conduit_pressure("steady_state_flow", iterations=100, file_prefix="constant_source_term")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

file_prefix = "constant_source_term"

plot_characteristics(folder, file_prefix)

We expected $\frac{dq}{dt} = S$ along the characteristic $\frac{dx}{dt} = u$. Let's verify that it true when the source is linear $S = 0.05$

# Let us plot a cross section along the charectaeristic that crosses the x-axis at x = -800m
def get_slope_along_characteristic(folder, file_prefix, u=2, N_m = 200, N_t = 100, T_end = 40, x_intercept = 200, x_length=1000):

    X = np.zeros((N_t, N_m)) 
    F = np.zeros((N_t, N_m)) 

    for i in range(0, N_t):
        solver = readwritedatafiles.read_data_file(f"{folder}/{file_prefix}_{i}.pkl")
        arhoX = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityX")]
        f, _, = physics.get_conv_flux_interior(solver.state_coeffs)

        X[i, :] = arhoX.ravel()
        F[i, :] = f[:,0,8,0].ravel()

    # define t and x values in steps
    t0 = 0
    t1 = N_t-1

    x0 = (x_intercept/x_length) * N_m
    x1 =((x_intercept + u* T_end)/x_length) * N_m 

    dq = (X[t1, int(x1)] - X[t0, int(x0)])
    dt = (T_end - 0)

    return dq/dt

# Let us plot a cross section along the charectaeristic that crosses the x-axis at x = -800m
def plot_dqdt_along_characteristic(folder, file_prefix, u=2, N_m = 200, N_t = 100, T_end = 40, x_intercept = 200, x_length=1000):

    X = np.zeros((N_t, N_m)) 
    F = np.zeros((N_t, N_m)) 

    dt = np.zeros(N_t)
    dq = np.zeros(N_t)

    for i in range(0, N_t):
        solver = readwritedatafiles.read_data_file(f"{folder}/{file_prefix}_{i}.pkl")
        arhoX = solver.state_coeffs[:,:,solver.physics.get_state_index("pDensityX")]
        f, _, = physics.get_conv_flux_interior(solver.state_coeffs)

        X[i, :] = arhoX.ravel()
        F[i, :] = f[:,0,8,0].ravel()

        dt[i] = i/N_t * T_end 
        # dx = u * dt 
        dq[i] = arhoX.ravel()[int((x_intercept + u * T_end)/x_length * N_m)]

    fig, ax = plt.subplots(figsize=(6,6))
    ax.plot(dt, dq)
    ax.set_xlabel("Time (s)")
    ax.set_ylabel("Q (new state variable)")
    ax.set_title("dq/dt along the characteristic that crosses the x-axis at x = -800m")

    plt.plot()
    

print(f"dq/dt = {get_slope_along_characteristic(folder, file_prefix)} along the characteristic that crosses the x-axis at x = -800m")

plot_dqdt_along_characteristic(folder, file_prefix)

dq/dt = 0.049553120069725246 along the characteristic that crosses the x-axis at x = -800m

4.2 Source term linear with time

I actually ran into some challenges when implementing this. At first glance the code mechanics of adding time dependence to the source term seemed quite simple. However, looking more closely it appears that time dependent source terms may not be supported.

  File "/Users/paxton/git/quail_volcano/src/solver/tools.py", line 318, in calculate_artificial_viscosity_integral
    Sq = physics.eval_source_terms(Uq, elem_helpers.x_elems,
        lambda: NotImplementedError(
                "Time dependence not implemented in Artificial Viscosity mod in tools.py"),
        Sq)
  File "/Users/paxton/git/quail_volcano/src/physics/base/base.py", line 680, in eval_source_terms
    Sq += source.get_source(self, Uq, xphys, time)
          ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/paxton/git/quail_volcano/src/physics/multiphasevpT/functions.py", line 4985, in get_source
    S[:, :, iarhoX] = 0.05 * t
                      ~~~~~^~~
TypeError: unsupported operand type(s) for *: 'float' and 'function'

I sent Fred an email and in the mean time I will try a source term dependent on position.

4.3 Source term linear with position

file_prefix = "linear_with_position_source_term"
ani = animate.animate_conduit_pressure("steady_state_flow", iterations=100, file_prefix=file_prefix)

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

plot_characteristics(folder, file_prefix)

4.4 Source term quadratic with position

file_prefix = "quadratic_with_position_source_term"
ani = animate.animate_conduit_pressure("steady_state_flow", iterations=100, file_prefix=file_prefix)

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test

plot_characteristics(folder, file_prefix)

5. Comparing solutions with different time step and space step size.

Let's compare the solution with sinusodial initial conditions with different time steps and step size.

Using 4000 time steps and 200 space steps we get the following.

file_prefix = "sinusoidal_IC_test"

plot_characteristics(folder, file_prefix)

If we increase our resolution by 4x for both space and time we now get this updated chart.

t_steps = 12,000

space_steps = 500

Higher spatial and temporal resolution

file_prefix = "sinusoidal_IC_test_high_res"

plot_characteristics(folder, file_prefix, N_m=500, N_t=100)

With the two charts above we can see that by increasing the spatial and temporal resolution we do have less diffusion. What happens if we just increase the temporal sampling?

Just higher temporal resolution:

file_prefix = "sinusoidal_IC_test_high_time_freq"

plot_characteristics(folder, file_prefix, N_m=200, N_t=100)

Alright that doesn't look so good. That implies that we are largely seeing wins because of smaller physical cells.

Just Higher spatial res

Ahh but it appears that they are coupled. In order to increase the spatial resolution to 500, I correspondingly had to increase the time resolution to basically 120,000 steps which is also what I plotted above.

6. Open questions

There were two things that I did not achieve with these experiments that I wanted to achieve.

1. I was unable to define time dependent source term because of the way the model is setup. I will followup with Frank on this one.
2. There is more diffusion in these simulations than I would like for the sake of clearly seeing the the initial conditions are correctly propagated along the characteristic lines. Perhaps I should look into overcoming the numerical stability issues to increase the solution order to 1?