solve_numpy.py

# translated from https://gitlab.com/manzhara157/nanobelts/-/blob/main/example.pas

import math
import os
import numpy
from contextlib import suppress


def clean_files():
    with suppress(FileNotFoundError):
        os.remove("mean_abc_numpy.txt")
    with suppress(FileNotFoundError):
        os.remove("fmean_v_numpy.txt")
    with suppress(FileNotFoundError):
        os.remove("fd_numpy.txt")
    with suppress(FileNotFoundError):
        os.remove("fn_numpy.txt")
    with suppress(FileNotFoundError):
        os.remove("fp_numpy.txt")
    with suppress(FileNotFoundError):
        os.remove("fv_numpy.txt")


def save_data(t, n, mean_a, mean_b, mean_c, mean_v, d, p, sumv):
    ln_t = math.log(t)
    with open("mean_abc_numpy.txt", "a") as f:
        f.write(
            "%f %f %f %f\n"
            % (ln_t, math.log(mean_a), math.log(mean_b), math.log(mean_c))
        )

    with open("fmean_v_numpy.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(mean_v)))

    with open("fd_numpy.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(d)))

    with open("fn_numpy.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(n)))

    with open("fp_numpy.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(p)))

    with open("fv_numpy.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(sumv)))


def solve():
    n0 = 200000

    z = 12

    tmax = 1000000
    dt = 0.1

    ra_kinet = 1
    rb_kinet = 1
    rc_kinet = 1

    ra_therm = 1
    rb_therm = 1
    rc_therm = 1

    ra_bal = 1
    rb_bal = 1
    rc_bal = 1

    j = 0.008
    xeq = 0.01
    d0 = 0.01

    a0 = 200.0
    b0 = 200.0
    c0 = 200.0
    vtot = n0 * a0 * b0 * c0 * 100

    lich0 = 100

    t = 0
    d = d0

    ra_bal = 1 / 2
    rb_bal = math.sqrt(2)
    rc_bal = math.sqrt(2)

    clean_files()

    # initialize data
    a = a0 * numpy.random.rand(n0)
    b = b0 * numpy.random.rand(n0)
    c = c0 * numpy.random.rand(n0)

    # a = numpy.fromfile("a.txt", sep="\n")
    # b = numpy.fromfile("b.txt", sep="\n")
    # c = numpy.fromfile("c.txt", sep="\n")

    # initial v

    v0 = numpy.sum(a * b * b)

    lich = lich0

    while t < tmax:
        t = t + dt

        anew = a + dt * ra_kinet * (
            d - ((rc_therm / c) + (rb_therm / b) + (ra_bal / ra_kinet) * j)
        )
        bnew = b + dt * rb_kinet * (
            d - ((ra_therm / a) + (rc_therm / c) + (rb_bal / rb_kinet) * j)
        )
        cnew = c + dt * rc_kinet * (
            d - ((ra_therm / a) + (rb_therm / b) + (rc_bal / rc_kinet) * j)
        )

        positives = numpy.logical_and(
            numpy.logical_and(anew >= 0, bnew >= 0), cnew >= 0
        )

        a = numpy.extract(positives, anew)
        b = numpy.extract(positives, bnew)
        c = numpy.extract(positives, cnew)

        sumv = numpy.sum(a * b * c)
        sump = numpy.sum(a * b + b * c + c * a)

        d = ((xeq + d0) * (1 - v0 / vtot) + v0 / vtot - sumv / vtot) / (
            1 - sumv / vtot
        ) - xeq
        p = 2 * sump

        if lich == lich0:
            n = a.size
            mean_a = numpy.mean(a)
            mean_b = numpy.mean(b)
            mean_c = numpy.mean(c)
            mean_v = sumv / n

            save_data(t, n, mean_a, mean_b, mean_c, mean_v, d, p, sumv)

            lich = 0

        lich = lich + 1


solve()