fipy_accuracy.py

from __future__ import division, print_function, absolute_import
import fipy as fp
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

axtickfsize = 18
labelfsize = 20
legfsize = labelfsize - 2
txtfsize = labelfsize - 2
lwidth = 3
markersize = 10
markeredgewidth = 0.1
mpl.rcParams['xtick.labelsize'] = axtickfsize
mpl.rcParams['ytick.labelsize'] = axtickfsize
mpl.rcParams['axes.labelsize'] = labelfsize
mpl.rcParams['font.size'] = txtfsize
mpl.rcParams['legend.fontsize'] = legfsize
mpl.rcParams['lines.linewidth'] = lwidth
mpl.rcParams['lines.markersize'] = markersize
mpl.rcParams['lines.markeredgewidth'] = markeredgewidth


def solve_and_find_error(dxvec, test_run=False):
    fpnp = fp.numerix
    mesh = fp.Grid1D(dx=dxvec)
    c = fp.CellVariable(name="concentration", mesh=mesh, value=0.)
    c_analyt = fp.CellVariable(name="analytical solution", mesh=mesh, value=0.)
    cntr = 0.5
    stretch = 1e-2
    C1 = 1/(fpnp.arcsinh(cntr/stretch) - fpnp.arcsinh(-cntr/stretch))
    C2 = -C1*fpnp.arcsinh(-cntr/stretch)
    c.constrain(0., mesh.facesLeft)
    c.constrain(1., mesh.facesRight)
#    print(mesh.cellVolumes)
#    print(dir(mesh))
#    zz
    X = mesh.x
    D = (1 + ((X - cntr)/stretch)**2)**(0.5)
    c_analyt = C1*fpnp.arcsinh((X - cntr)/stretch) + C2
    eq = 0 == fp.DiffusionTerm(coeff=D, var=c)
    eq.solve()
#    error = (fpnp.sum((c - c_analyt)**2)/len(c))**(0.5)
    error = (fpnp.sum((c - c_analyt)**2 * mesh.cellVolumes))**(0.5)
    if test_run:
        viewer = fp.Viewer(vars=(c, c_analyt), datamin=0., datamax=1.)
        viewer.plot()
        input()
    return error


def dxvec_nonuniform_exp(nx):
    reduce_factor = 0.98
    dx0 = 1/(sum([reduce_factor**i for i in range(nx)]))
    dxvec = np.array([dx0*reduce_factor**i for i in range(nx)])
    return dxvec


def dxvec_nonuniform_segments(nx):
    L1, L2, L3 = 0.4, 0.2, 0.4
    n1, n3 = int(0.15*nx), int(0.1*nx)
    n2 = nx - n1 - n3
    dx1, dx2, dx3 = L1/n1, L2/n2, L3/n3
    dxvec = [dx1]*n1 + [dx2]*n2 + [dx3]*n3
    return dxvec


def main():
    nxvec = np.arange(150, 200)
    errors_u = np.zeros(len(nxvec))
    errors_nu_seg = np.zeros(len(nxvec))
    errors_nu_exp = np.zeros(len(nxvec))
#    # TEST
#    nx = 700
#    dxvec = dxvec_nonuniform_exp(nx)
#    solve_and_find_error(dxvec, test_run=True)
#    #
    for indx, nx in enumerate(nxvec):
        dxvec_u = [1/nx]*nx
        dxvec_nu_seg = dxvec_nonuniform_segments(nx)
        dxvec_nu_exp = dxvec_nonuniform_exp(nx)
        errors_u[indx] = solve_and_find_error(dxvec_u)
        errors_nu_seg[indx] = solve_and_find_error(dxvec_nu_seg)
        errors_nu_exp[indx] = solve_and_find_error(dxvec_nu_exp)
    log_nxvec = np.log(nxvec)
    log_errors_u = np.log(errors_u)
    log_errors_nu_seg = np.log(errors_nu_seg)
    log_errors_nu_exp = np.log(errors_nu_exp)
    fit_u = np.polyfit(log_nxvec, log_errors_u, 1)
    fit_nu_seg = np.polyfit(log_nxvec, log_errors_nu_seg, 1)
    fit_nu_exp = np.polyfit(log_nxvec, log_errors_nu_exp, 1)

    fig, axs = plt.subplots(1, 3, figsize=(16, 9))
    ax = axs[0]
    ax.plot(log_nxvec, log_errors_u, label="simulation errors")
    ax.plot(log_nxvec, np.poly1d(fit_u)(log_nxvec), linestyle="--",
            label="Line: {m:1.2f}x+{b:1.2f}".format(m=fit_u[0], b=fit_u[1]))
    ax.legend(loc="best")
    ax.set_xlabel("log(nx)")
    ax.set_ylabel("log(error)")
    ax.set_title("Uniform mesh")

    ax = axs[1]
    ax.plot(log_nxvec, log_errors_nu_seg, label="simulation errors")
    ax.plot(log_nxvec, np.poly1d(fit_nu_seg)(log_nxvec), linestyle="--",
            label="Line: {m:1.2f}x+{b:1.2f}".format(m=fit_nu_seg[0], b=fit_nu_seg[1]))
    ax.legend(loc="best")
    ax.set_xlabel("log(nx)")
    ax.set_ylabel("log(error)")
    ax.set_title("Non-uniform (segments) mesh")

    ax = axs[2]
    ax.plot(log_nxvec, log_errors_nu_exp, label="simulation errors")
    ax.plot(log_nxvec, np.poly1d(fit_nu_exp)(log_nxvec), linestyle="--",
            label="Line: {m:1.2f}x+{b:1.2f}".format(m=fit_nu_exp[0], b=fit_nu_exp[1]))
    ax.legend(loc="best")
    ax.set_xlabel("log(nx)")
    ax.set_ylabel("log(error)")
    ax.set_title("Non-uniform (exp) mesh")

    fig.savefig("convergence.pdf", bbox_inches="tight")
    plt.show()
    return

if __name__ == "__main__":
    main()