#!/usr/bin/env python3
"``naginterfaces.library.pde.dim2_gen_order2_rectangle`` Python Example."
# NAG Copyright 2020.
# pylint: disable=invalid-name,too-many-locals,too-many-arguments
from math import pi, sin
import numpy as np
from naginterfaces.library import machine, pde
[docs]def main():
"""
Example for :func:`naginterfaces.library.pde.dim2_gen_order2_rectangle`.
General system of second-order PDEs, method of lines, finite differences,
remeshing, two space variables, rectangular region.
>>> main()
naginterfaces.library.pde.dim2_gen_order2_rectangle Python Example Results.
Solve a predator-prey problem.
Solution at every 2nd grid point in level 3 at time 0.0250:
x y approx c1 approx c2
...
1.000e+00 1.000e+00 1.132e+02 1.132e+06
"""
print(
'naginterfaces.library.pde.dim2_gen_order2_rectangle '
'Python Example Results.'
)
print('Solve a predator-prey problem.')
ts = 0.
twant = [0.01, 0.025]
dt = [0.5e-3, 1.e-7, 0.]
xmin = 0.
xmax = 1.
ymin = 0.
ymax = 1.
nx = 11
ny = 11
tols = 0.1
tolt = 0.1
alpha = 50.
beta = 300.
def cb_pdedef(t, x, y, u, ut, ux, uy, uxx, uxy, uyy):
"""The system of PDEs."""
res = np.empty(u.shape)
fourpi = 4.*pi
for i, xi in enumerate(x):
b1 = (
1. + alpha*xi*y[i] +
beta*sin(fourpi*xi)*sin(fourpi*y[i])
)
res[i, 0] = (
ut[i, 0] - (uxx[i, 0]+uyy[i, 0]) -
u[i, 0]*(b1-u[i, 0]-0.5E-6*u[i, 1])
)
res[i, 1] = (
-0.05*(uxx[i, 1]+uyy[i, 1]) -
u[i, 1]*(-b1+1.0E4*u[i, 0]-u[i, 1])
)
return res
def cb_bndary(t, x, y, u, ut, ux, uy, lbnd, res):
"""The boundary conditions."""
tol = 10.*machine.precision()
new_res = np.empty(u.shape)
new_res[:, :] = res
for jp in lbnd:
j = jp - 1
if abs(x[j]) <= tol or abs(x[j]-1.) <= tol:
new_res[j, :] = ux[j, :]
elif abs(y[j]) <= tol or abs(y[j]-1.) <= tol:
new_res[j, :] = uy[j, :]
return new_res
def cb_pdeiv(npde, t, x, y):
"""The PDE initial values."""
u = np.empty((len(x), npde))
fourpi = 4.*pi
u[:, 0] = 10. + (16.*x*(1.-x)*y*(1.-y))**2
u[:, 1] = (
-1. - alpha*x*y -
beta*np.sin(fourpi*x)*np.sin(fourpi*y) + 1.0E4*u[:, 0]
)
return u
def cb_monitr_dummy(
npde, t, dt, dtnew, tlast, ngpts, xpts, ypts, lsol, sol, ierr,
):
"""Dummy monitoring function."""
return ierr
def cb_monitr(
npde, t, dt, dtnew, tlast, ngpts, xpts, ypts, lsol, sol, ierr,
):
"""Monitoring function."""
if not tlast:
return ierr
nlev = len(ngpts)
level = nlev
print(
'Solution at every 2nd grid point ' +
'in level {:d} at time {:8.4f}:'.format(level, t)
)
npts = ngpts[level-1]
ipsol = lsol[level-1]
k = sum(ngpts[:nlev-1])
print(
' '*7 + 'x' + ' '*10 + 'y' +
' '*9 + 'approx c1' + ' '*4 + 'approx c2'
)
for i in range(0, npts, 2):
print(' {:11.3e} {:11.3e} {:11.3e} {:11.3e}'.format(
xpts[k+i], ypts[k+i], sol[ipsol+i], sol[ipsol+npts+i],
))
return ierr
opti = [4, 0, 0, 0]
optr = [
[250., 1.5e6],
[1., 1.],
[1., 1.],
]
comm = {}
itrace = 0
ind = 0
for t_i, tout in enumerate(twant):
ts, dt, ind = pde.dim2_gen_order2_rectangle(
ts, tout, dt, xmin, xmax, ymin, ymax, nx, ny, tols, tolt,
cb_pdedef, cb_bndary, cb_pdeiv,
cb_monitr_dummy if t_i == 0 else cb_monitr,
opti, optr, comm, itrace, ind,
)
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(
None, verbose=True, report=False,
optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
).failed
)
