Last active
January 18, 2018 17:28
-
-
Save raybsmith/e57f6f4739e24ff9c97039ad573a3621 to your computer and use it in GitHub Desktop.
Test accuracy of fipy with non-uniform mesh
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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): | |
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) | |
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) | |
return error | |
def dxvec_nonuniform_exp(nx): | |
reduce_factor = 0.96 | |
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)) | |
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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment