from ChemPlugin import *
import sys

def write_line(f, cp, gap, then):
    now = cp[0].Report1("Time", "years")
    if (then - now) <= gap / 1e4:
        # Writing Temperature profile to file
        f.write("%8.4f" % now)
        for c in cp:
            f.write("\t%8.4f" % c.Report1("temperature", "C"))
        f.write("\n")
        then += gap
    return then

def end_run(code):
    input()
    sys.exit(code)

print("Model heat conduction in one dimension\n")

# Simulation Parameters
nx = 100                     # number of instances along x
length = 100.0               # m
deltax = length / nx         # m
deltay = 1.0                 # m
deltaz = 1.0                 # m
porosity = 0.25              # volume fraction

tcond = 2.0                  # W/m/k

time_end = 15.0              # years
delta_years = time_end / 3.0 # years
next_output = 0.0            # years

# Open output file and write instance positions on the first line.
f = open("HeatConduction.txt", 'w')
f.write("years")
for i in range(nx):
    f.write("\t%8.4f" % ((i + 0.5) * deltax))
f.write("\n")

# Configure and Initialize the instances
cp = [ChemPlugin() for i in range(nx)]

cp[0].Console("stdout")
cp[0].Config("pluses = banner")

cmd = "span 20 C to 100 C; " +\
      "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg; " +\
      "volume = " + str(deltax * deltay * deltaz) + " m3;" +\
      "porosity = " + str(porosity) + ";" +\
      "time end = " + str(time_end) + " years"

for i in range(nx):
    cp[i].Config(cmd)
    if (i < nx / 2):
        cp[i].Config("T = 100 C")
    else:
        cp[i].Config("T = 20 C")

    if cp[i].Initialize():
        print("Instance " + str(i) + " failed to initialize")
        end_run(-1)

next_output = write_line(f, cp, delta_years, next_output)

# Link the instances
trans = deltay * deltaz * tcond / deltax  # W/K
for i in range(1, nx):
    link = cp[i].Link(cp[i - 1])
    link.HeatTrans(trans, "W/K")

while True:
    deltat = 1e99
    for c in cp:
        deltat = min(deltat, c.ReportTimeStep())
    for c in cp:
        if c.AdvanceTimeStep(deltat): end_run(0)
    for c in cp:
        if c.AdvanceHeatTransport(): end_run(-1)
    for c in cp:
        if c.AdvanceChemical(): end_run(-1)

    next_output = write_line(f, cp, delta_years, next_output)