Phase Appearance/Disappearance

import os
from pathlib import Path

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

import numpy as np
from ogs6py.ogs import OGS

model = OGS(INPUT_FILE="bourgeat.prj", PROJECT_FILE=f"{out_dir}/modified.prj")
# This Jupyter notebook version of this test runs not as far as its cTest-counterpart,
# it'll stop after approx. 800 ka
model.replace_text(2.5e13, xpath="./time_loop/processes/process/time_stepping/t_end")


# The cTest version shows only a few output-timesteps while this version is supposed to
# output much higher resolution time steps to be able to compare the results. Thus, every
# timestep will be written and the maximum timestep-size of the adaptive time stepping
# method will be reduced drastically
time_end = 1e11
model.replace_text(
    time_end, xpath="./time_loop/processes/process/time_stepping/maximum_dt"
)

# The following for loop generates a text with output times,which is then replaced by
# the ogs6py API in the project file.
new_line = "\u000A"
timesteps = str(0.4 * 1e6 * 86400 * 365.25) + new_line
for t in np.arange(0.6, 1.1, 0.1):
    timesteps += str(t * 1e6 * 86400 * 365.25) + new_line

model.replace_text(1, xpath="./time_loop/output/timesteps/pair/each_steps")
model.replace_text(timesteps, xpath="./time_loop/output/fixed_output_times")
model.write_input()

# Run OGS
model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")

# Colors
cls = ["#e6191d", "#337fb8", "#4eae4c", "#984ea3", "#984ea3", "#feff32"]

import vtuIO

# Read PVD-output
pvdfile = vtuIO.PVDIO(f"{out_dir}/result_bourgeat.pvd", dim=2)
point = {"A": (0.0, 0.0, 0.0)}
time = pvdfile.timesteps

saturation = pvdfile.read_time_series("saturation", point)
gas_pressure = pvdfile.read_time_series("gas_pressure", point)
liquid_pressure = pvdfile.read_time_series("liquid_pressure_interpolated", point)

num_results = [1.0 - saturation["A"], gas_pressure["A"], liquid_pressure["A"]]

time_years = time / 365.2425 / 86400

import pandas as pd

# Read the reference data from CSV files
refs = [
    pd.read_csv("references/bourgeat_sG.csv"),
    pd.read_csv("references/bourgeat_pGR.csv"),
    pd.read_csv("references/bourgeat_pLR.csv"),
]

header = list(refs[0].keys())

indices = {"Gas saturation": 0, "Gas pressure": 1, "Liquid pressure": 2}
labels = ["$s_{G}$", "$p_{GR}$", "$p_{LR}$"]

import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = (12, 4)

# Loop over gas_saturation, gas_pressure, and liquid_pressure
for i in indices:
    index = indices[i]
    ref_time = refs[index]["time"]

    fig1, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
    fig1.suptitle(i + r" vs. time at $x=0$m")

    ax1.set_xscale("log")
    ax1.set_xlabel("time / a", fontsize=12)
    ax1.set_ylabel(labels[index], fontsize=12)

    for r in range(1, len(refs[index].columns)):
        ax1.plot(
            ref_time,
            refs[index][refs[index].keys()[r]],
            linewidth=1,
            linestyle="-",
            label=refs[index].keys()[r],
        )

    ax1.plot(
        time_years,
        num_results[index],
        color="black",
        linewidth=2,
        linestyle="--",
        label="OGS-TH$^2$M",
    )

    ax2.set_xlabel("time / a", fontsize=12)

    for r in range(1, len(refs[index].columns)):
        ax2.plot(
            ref_time,
            refs[index][refs[index].keys()[r]],
            linewidth=1,
            linestyle="-",
            label=refs[index].keys()[r],
        )
    ax2.plot(
        time_years,
        num_results[index],
        color="black",
        linewidth=2,
        linestyle="--",
        label="OGS-TH$^2$M",
    )

    ax1.legend()


fig1.savefig("results_sG_pGR_pLR.pdf")