Source code for naginterfaces.library.examples.pde.dim2_gen_order2_rectangle_ex

#!/usr/bin/env python
"``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 )