02_create_run_prj.py

"""
Summer-school example (2/3): create the OGS project file and run it.

Builds the OpenGeoSys `HEAT_TRANSPORT_BHE` project (.prj) for the 37-BHE BTES
mesh produced by `01_create_mesh.py`, using the ogstools Project (ogs6py) API,
then runs OGS.

The mesh (./model/) provides the MaterialIDs the project file maps onto:
  * MaterialID 0 -> weathered basalt (soil medium 0)
  * MaterialID 1 -> granodiorite     (soil medium 1)
  * MaterialIDs 2..38 -> the 37 BHEs (one borehole_heat_exchanger + one
    temperature_BHE process variable each)

The BHEs are coaxial (CXA) and run on a seasonal cycle (inlet-temperature
control): inject/charge at 90 degC for the first half year, then extract/
discharge at 30 degC for the second half year, at a constant 4 L/s per BHE.
This seasonal cycle is repeated for 5 years; the adaptive time stepper is forced
to land exactly on every inject<->extract switch (fixed_output_times).
The initial soil temperature follows a geothermal gradient, held by Dirichlet
temperatures on the top (surface) and bottom faces; the sides are insulated.

NOTE: all material/operational values below are SENSIBLE PLACEHOLDERS for a
workshop demo - replace them with your own data for a real study.
"""

import subprocess
import sys
from pathlib import Path

import ogstools as ot

# ===========================================================================
# PARAMETERS  --  edit these
# ===========================================================================
HERE = Path(__file__).resolve().parent
MODEL_DIR = HERE / "model"                 # holds the meshes + the .prj
RESULTS_DIR = MODEL_DIR / "results"        # OGS writes output here
PRJ_FILE = MODEL_DIR / "btes_model.prj"
RUN_OGS = True                             # False -> only write the .prj

N_BHE = 37
# Insulated sides are the natural zero-flux BC in OGS, so no 'sides' mesh.
MESHES = ["domain", "top", "bottom"]

# --- Geometry (MUST match 01_create_mesh.py) -------------------------------
DEPTH = 500.0            # total domain depth [m]
BHE_LENGTH = 400.0       # BHE length [m]  (== bhe_length in the mesh)
BHE_DIAMETER = 0.2413    # borehole diameter [m]  (== 2 * bhe_radius)

# --- Geothermal field (degC): initial + boundary temperatures --------------
SURFACE_TEMP = 10.0                                # ground-surface temperature
GEO_GRADIENT = 0.03                                # geothermal gradient [K/m]
BOTTOM_TEMP = SURFACE_TEMP + GEO_GRADIENT * DEPTH  # temperature at z = -DEPTH
BHE_INIT_TEMP = SURFACE_TEMP + GEO_GRADIENT * BHE_LENGTH / 2.0

# --- Coaxial (CXA) BHE: seasonal operation (inlet-temperature control) ------
BHE_TYPE = "CXA"          # "CXA" (annular inlet) or "CXC" (centred inlet)
BHE_FLOW_RATE = 0.004     # [m3/s] = 4 L/s per BHE, constant for the whole run
INJECT_TEMP = 90.0        # [degC] inlet temperature while injecting (1st half year)
EXTRACT_TEMP = 30.0       # [degC] inlet temperature while extracting (2nd half year)

# --- Coaxial pipe geometry (PLACEHOLDER) -----------------------------------
# outer = pipe at the borehole wall, inner = central pipe
PIPE_OUTER = dict(diameter=0.18, wall_thickness=0.006, wall_k=1.3)
PIPE_INNER = dict(diameter=0.10, wall_thickness=0.007, wall_k=0.4)
LONG_DISP_LENGTH = 0.001  # longitudinal dispersion length [m]

# --- Grout, refrigerant (PLACEHOLDER) --------------------------------------
GROUT = dict(density=2190.0, porosity=0.1, cp=735.16, k=2.2)
REFRIGERANT = dict(density=1000.0, viscosity=0.00114, cp=4190.0, k=0.59,
                   reference_temperature=25.0)

# --- Soil media (PLACEHOLDER thermal properties) ---------------------------
# Convention (as in the OGS BHE examples): the Solid phase carries the
# volumetric heat capacity via density=1 and specific_heat_capacity = rho*cp.
# MaterialID 0 = weathered basalt (0..-200 m), MaterialID 1 = granodiorite.
WEATHERED_BASALT = dict(porosity=0.15, thermal_conductivity=1.7,
                        solid_vol_heat_capacity=2.3e6, darcy_velocity="0 0 0")
GRANODIORITE = dict(porosity=0.01, thermal_conductivity=3.1,
                    solid_vol_heat_capacity=2.4e6, darcy_velocity="0 0 0")

# --- Operation period & time stepping --------------------------------------
N_YEARS = 5                   # number of seasonal cycles to simulate
YEAR = 365 * 86400.0          # [s] one seasonal cycle = 1 year
T_END = N_YEARS * YEAR        # [s] total simulated time (5 years)
SWITCH_TIME = YEAR / 2.0      # [s] inject -> extract switch within each year
INLET_RAMP = 86400.0          # [s] ramp window for each inlet-temperature switch

# Adaptive (iteration-number-based) time stepping: dt grows when the nonlinear
# solver converges easily and shrinks at the stiff transients (injection start
# and every seasonal switch). A step that fails to converge is rejected and
# retried with a smaller dt automatically.
INITIAL_DT = 600.0           # [s] first step (10 min)
MIN_DT = 60.0                # [s] lower bound (1 min)
MAX_DT = 5 * 86400.0         # [s] upper bound (5 days)
# Paired entrywise: if the last step took <= ADAPT_ITERS[k] nonlinear iterations,
# multiply dt by ADAPT_MULT[k] for the next step.
ADAPT_ITERS = [6, 10, 20, 40]
ADAPT_MULT = [2.0, 1.6, 0.8, 0.4]

OUTPUT_PREFIX = "btes"         # output every computed time step (each dt) to VTU
# ===========================================================================

bhe_pvs = [f"temperature_BHE{i}" for i in range(1, N_BHE + 1)]

prj = ot.Project(output_file=PRJ_FILE)

# --- meshes ----------------------------------------------------------------
for name in MESHES:
    prj.mesh.add_mesh(filename=f"{name}.vtu")

# --- process + process variables -------------------------------------------
prj.processes.set_process(
    name="HeatTransportBHE", type="HEAT_TRANSPORT_BHE", integration_order="2"
)
prj.processes.add_process_variable(
    process_variable="process_variable", process_variable_name="temperature_soil"
)
for name in bhe_pvs:
    prj.processes.add_process_variable(
        process_variable="process_variable", process_variable_name=name
    )

# --- borehole heat exchangers ----------------------------------------------
# ogstools' high-level BHE builder targets an older OGS schema (it emits the
# removed "FixedPowerConstantFlow" control and a coaxial "distance_between_pipes"
# that OGS 6.5.x rejects), so we write the current form with the generic
# add_element/add_block methods. All 37 BHEs are identical, so a single
# add_* call with the borehole_heat_exchanger xpath fills every one of them.
proc = "./processes/process"
prj.add_element(parent_xpath=proc, tag="borehole_heat_exchangers")
bhes = f"{proc}/borehole_heat_exchangers"
for i in range(N_BHE):
    prj.add_element(parent_xpath=bhes, tag="borehole_heat_exchanger",
                    attrib_list=["id"], attrib_value_list=[str(i)])

bhe = f"{bhes}/borehole_heat_exchanger"          # matches all 37
prj.add_element(parent_xpath=bhe, tag="type", text=BHE_TYPE)
# inlet-temperature control: prescribe the inlet temperature (via the time-curve
# parameter BHE_inflow_temp) together with a constant flow rate.
prj.add_block("flow_and_temperature_control", parent_xpath=bhe,
              taglist=["type", "flow_rate", "temperature"],
              textlist=["InflowTemperature", BHE_FLOW_RATE, "BHE_inflow_temp"])
prj.add_block("borehole", parent_xpath=bhe,
              taglist=["length", "diameter"], textlist=[BHE_LENGTH, BHE_DIAMETER])
prj.add_block("grout", parent_xpath=bhe,
              taglist=["density", "porosity", "specific_heat_capacity", "thermal_conductivity"],
              textlist=[GROUT["density"], GROUT["porosity"], GROUT["cp"], GROUT["k"]])
# coaxial pipes: <outer>/<inner> (no distance_between_pipes)
prj.add_element(parent_xpath=bhe, tag="pipes")
pipes = f"{bhe}/pipes"
prj.add_block("outer", parent_xpath=pipes,
              taglist=["diameter", "wall_thickness", "wall_thermal_conductivity"],
              textlist=[PIPE_OUTER["diameter"], PIPE_OUTER["wall_thickness"], PIPE_OUTER["wall_k"]])
prj.add_block("inner", parent_xpath=pipes,
              taglist=["diameter", "wall_thickness", "wall_thermal_conductivity"],
              textlist=[PIPE_INNER["diameter"], PIPE_INNER["wall_thickness"], PIPE_INNER["wall_k"]])
prj.add_element(parent_xpath=pipes, tag="longitudinal_dispersion_length", text=LONG_DISP_LENGTH)
prj.add_block("refrigerant", parent_xpath=bhe,
              taglist=["density", "viscosity", "specific_heat_capacity",
                       "thermal_conductivity", "reference_temperature"],
              textlist=[REFRIGERANT["density"], REFRIGERANT["viscosity"], REFRIGERANT["cp"],
                        REFRIGERANT["k"], REFRIGERANT["reference_temperature"]])

# --- media (one per soil MaterialID; BHE materials need no medium) ----------
for mid, props in ((0, WEATHERED_BASALT), (1, GRANODIORITE)):
    prj.media.add_property(medium_id=mid, phase_type="Solid",
                           name="specific_heat_capacity", type="Constant",
                           value=props["solid_vol_heat_capacity"])
    prj.media.add_property(medium_id=mid, phase_type="Solid",
                           name="density", type="Constant", value=1)
    prj.media.add_property(medium_id=mid, phase_type="AqueousLiquid",
                           name="phase_velocity", type="Constant",
                           value=props["darcy_velocity"])
    prj.media.add_property(medium_id=mid, phase_type="AqueousLiquid",
                           name="specific_heat_capacity", type="Constant", value=4190)
    prj.media.add_property(medium_id=mid, phase_type="AqueousLiquid",
                           name="density", type="Constant", value=1000)
    prj.media.add_property(medium_id=mid, name="porosity", type="Constant",
                           value=props["porosity"])
    prj.media.add_property(medium_id=mid, name="thermal_conductivity",
                           type="Constant", value=props["thermal_conductivity"])
    prj.media.add_property(medium_id=mid, name="thermal_longitudinal_dispersivity",
                           type="Constant", value=0)
    prj.media.add_property(medium_id=mid, name="thermal_transversal_dispersivity",
                           type="Constant", value=0)

# --- parameters ------------------------------------------------------------
# Geothermal-gradient initial soil temperature: T(z) = T_surface - gradient*z
# (z is negative downward, so temperature increases with depth).
prj.parameters.add_parameter(name="T0_soil", type="Function",
                             expression=f"{SURFACE_TEMP} - {GEO_GRADIENT}*z")
prj.parameters.add_parameter(name="T_surface", type="Constant", value=SURFACE_TEMP)
prj.parameters.add_parameter(name="T_bottom", type="Constant", value=BOTTOM_TEMP)
prj.parameters.add_parameter(name="T0_BHE", type="Constant",
                             values=" ".join([str(BHE_INIT_TEMP)] * 3))
# Seasonal BHE inlet temperature, repeated for N_YEARS: 90 degC (inject) for the
# first half of each year, then 30 degC (extract) for the second half, with a
# short INLET_RAMP at every switch so the inflow temperature stays continuous.
# The switch times are collected in `switch_times` so the time stepper can be
# forced to land exactly on them (fixed_output_times, below). InflowTemperature
# reads the curve through the CurveScaled parameter BHE_inflow_temp.
inflow_coords = [0.0]
inflow_values = [INJECT_TEMP]
switch_times = []
for k in range(N_YEARS):
    y0 = k * YEAR
    t_down = y0 + SWITCH_TIME                 # inject -> extract (mid-year)
    inflow_coords += [t_down, t_down + INLET_RAMP]
    inflow_values += [INJECT_TEMP, EXTRACT_TEMP]
    switch_times += [t_down, t_down + INLET_RAMP]
    t_up = y0 + YEAR                          # extract -> inject (year boundary)
    if k < N_YEARS - 1:
        inflow_coords += [t_up, t_up + INLET_RAMP]
        inflow_values += [EXTRACT_TEMP, INJECT_TEMP]
        switch_times += [t_up, t_up + INLET_RAMP]
    else:
        inflow_coords += [t_up]               # final point closes the curve at T_END
        inflow_values += [EXTRACT_TEMP]

prj.curves.add_curve(name="inflow_temperature",
                     coords=inflow_coords, values=inflow_values)
prj.parameters.add_parameter(name="Curve_Scale", type="Constant", value=1)
prj.parameters.add_parameter(name="BHE_inflow_temp", type="CurveScaled",
                             curve="inflow_temperature", parameter="Curve_Scale")

# --- process-variable definitions (ICs + BCs) ------------------------------
prj.process_variables.set_ic(process_variable_name="temperature_soil",
                       components="1", order="1", initial_condition="T0_soil")
prj.process_variables.add_bc(process_variable_name="temperature_soil", type="Dirichlet",
                       mesh="top", parameter="T_surface")
prj.process_variables.add_bc(process_variable_name="temperature_soil", type="Dirichlet",
                       mesh="bottom", parameter="T_bottom")
for name in bhe_pvs:
    prj.process_variables.set_ic(process_variable_name=name, components="3",
                           order="1", initial_condition="T0_BHE")

# --- time loop -------------------------------------------------------------
# Output every computed time step (each dt). In addition, every seasonal switch
# time is passed as fixed_output_times: this forces the adaptive stepper to
# shorten dt as needed so it lands exactly on each inject<->extract switch,
# resolving the seasonal forcing sharply instead of stepping over it.

prj.time_loop.add_process(process="HeatTransportBHE",
                          nonlinear_solver_name="basic_picard",
                          convergence_type="DeltaX", norm_type="NORM2",
                          reltol="1e-3", time_discretization="BackwardEuler")
prj.time_loop.set_stepping(
    process="HeatTransportBHE", type="IterationNumberBasedTimeStepping",
    t_initial="0", t_end=str(int(T_END)),
    initial_dt=str(int(INITIAL_DT)), minimum_dt=str(int(MIN_DT)),
    maximum_dt=str(int(MAX_DT)),
    number_iterations=ADAPT_ITERS, multiplier=ADAPT_MULT,
)
prj.time_loop.add_output(type="VTK", prefix=OUTPUT_PREFIX,
                         suffix="_ts_{:timestep}_t_{:time}",
                         repeat="1000000", each_steps="1",
                         fixed_output_times=switch_times,
                         variables=["temperature_soil"] + bhe_pvs)

# --- solvers ---------------------------------------------------------------
prj.nonlinear_solvers.add_non_lin_solver(name="basic_picard", type="Picard",
                                         max_iter="100", linear_solver="general_linear_solver")
prj.linear_solvers.add_lin_solver(name="general_linear_solver", kind="eigen",
                                  solver_type="BiCGSTAB", precon_type="ILUT",
                                  max_iteration_step="10000", error_tolerance="1e-10")

# --- write -----------------------------------------------------------------
prj.write_input()
print(f"Wrote project file: {PRJ_FILE}")
print(f"  {N_BHE} {BHE_TYPE} BHEs, length {BHE_LENGTH} m, diameter {BHE_DIAMETER} m")
print(f"  geothermal IC: {SURFACE_TEMP} degC -> {BOTTOM_TEMP:.0f} degC at -{DEPTH:.0f} m")
print(f"  operation: {N_YEARS} yearly cycles, inject {INJECT_TEMP:.0f} degC (half year) "
      f"then extract {EXTRACT_TEMP:.0f} degC (half year), flow {BHE_FLOW_RATE * 1000:.0f} L/s/BHE")
print(f"  time: adaptive dt (start {INITIAL_DT / 3600:.1f} h, max {MAX_DT / 86400:.0f} d), "
      f"t_end = {T_END:.0f} s ({N_YEARS} yr); output every step, "
      f"landing on {len(switch_times)} switch times")

# --- run OGS (live output in this terminal) --------------------------------
# Note: prj.run_model() either writes stdout to a log file (write_logs=True) or
# discards it to DEVNULL (write_logs=False) - it cannot stream to the terminal.
# So we call the OGS executable directly and let it inherit stdout/stderr, which
# shows OGS's progress live as it runs.
if RUN_OGS:
    RESULTS_DIR.mkdir(parents=True, exist_ok=True)
    ogs_exe = Path(sys.executable).parent / ("ogs.exe" if sys.platform == "win32" else "ogs")
    cmd = [str(ogs_exe), str(PRJ_FILE), "-o", str(RESULTS_DIR)]
    print("Running:", " ".join(cmd), flush=True)
    subprocess.run(cmd, check=True)     # inherits the terminal -> live output
    print("OGS run finished. Output written to", RESULTS_DIR)