Source code for pyro.multigrid.examples.mg_test_general_alphabeta_only

#!/usr/bin/env python3

"""Test the general MG solver with a variable coefficient Helmholtz
problem.  This ensures we didn't screw up the base functionality here.

Here we solve::

   alpha phi + div . ( beta grad phi ) = f

with::

   alpha = 1.0
   beta = 2.0 + cos(2.0*pi*x)*cos(2.0*pi*y)

   f = (-16.0*pi**2*cos(2*pi*x)*cos(2*pi*y) - 16.0*pi**2 + 1.0)*sin(2*pi*x)*sin(2*pi*y)


This has the exact solution::

   phi = sin(2.0*pi*x)*sin(2.0*pi*y)

on [0,1] x [0,1]

We use Dirichlet BCs on phi.  For beta, we do not have to impose the
same BCs, since that may represent a different physical quantity.
Here we take beta to have Neumann BCs.  (Dirichlet BCs for beta will
force it to 0 on the boundary, which is not correct here)

"""


import matplotlib.pyplot as plt
import numpy as np

import pyro.mesh.boundary as bnd
import pyro.multigrid.general_MG as MG
import pyro.util.io_pyro as io
from pyro.mesh import patch
from pyro.util import compare, msg


# the analytic solution
[docs] def true(x, y): return np.sin(2.0*np.pi*x)*np.sin(2.0*np.pi*y)
# the coefficients
[docs] def alpha(x, _y): return np.ones_like(x)
[docs] def beta(x, y): return 2.0 + np.cos(2.0*np.pi*x)*np.cos(2.0*np.pi*y)
[docs] def gamma_x(x, _y): return np.zeros_like(x)
[docs] def gamma_y(x, _y): return np.zeros_like(x)
# the righthand side
[docs] def f(x, y): return (-16.0*np.pi**2*np.cos(2*np.pi*x)*np.cos(2*np.pi*y) - 16.0*np.pi**2 + 1.0)*np.sin(2*np.pi*x)*np.sin(2*np.pi*y)
[docs] def test_general_poisson_dirichlet(N, store_bench=False, comp_bench=False, bench_dir="tests/", make_plot=False, verbose=1, rtol=1.e-12): """ test the general MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark """ # test the multigrid solver nx = N ny = nx # create the coefficient variable g = patch.Grid2d(nx, ny, ng=1) d = patch.CellCenterData2d(g) bc_c = bnd.BC(xlb="neumann", xrb="neumann", ylb="neumann", yrb="neumann") d.register_var("alpha", bc_c) d.register_var("beta", bc_c) d.register_var("gamma_x", bc_c) d.register_var("gamma_y", bc_c) d.create() a = d.get_var("alpha") a[:, :] = alpha(g.x2d, g.y2d) b = d.get_var("beta") b[:, :] = beta(g.x2d, g.y2d) gx = d.get_var("gamma_x") gx[:, :] = gamma_x(g.x2d, g.y2d) gy = d.get_var("gamma_y") gy[:, :] = gamma_y(g.x2d, g.y2d) # create the multigrid object a = MG.GeneralMG2d(nx, ny, xl_BC_type="dirichlet", yl_BC_type="dirichlet", xr_BC_type="dirichlet", yr_BC_type="dirichlet", coeffs=d, verbose=verbose, vis=0, true_function=true) # initialize the solution to 0 a.init_zeros() # initialize the RHS using the function f rhs = f(a.x2d, a.y2d) a.init_RHS(rhs) # solve to a relative tolerance of 1.e-11 a.solve(rtol=1.e-11) # alternately, we can just use smoothing by uncommenting the following # a.smooth(a.nlevels-1,50000) # get the solution v = a.get_solution() # compute the error from the analytic solution b = true(a.x2d, a.y2d) e = v - b enorm = e.norm() print(" L2 error from true solution = %g\n rel. err from previous cycle = %g\n num. cycles = %d" % (enorm, a.relative_error, a.num_cycles)) # plot the solution if make_plot: plt.clf() plt.figure(figsize=(10.0, 4.0), dpi=100, facecolor='w') plt.subplot(121) img1 = plt.imshow(np.transpose(v.v()), interpolation="nearest", origin="lower", extent=[a.xmin, a.xmax, a.ymin, a.ymax]) plt.xlabel("x") plt.ylabel("y") plt.title(f"nx = {nx}") plt.colorbar(img1) plt.subplot(122) img2 = plt.imshow(np.transpose(e.v()), interpolation="nearest", origin="lower", extent=[a.xmin, a.xmax, a.ymin, a.ymax]) plt.xlabel("x") plt.ylabel("y") plt.title("error") plt.colorbar(img2) plt.tight_layout() plt.savefig("mg_general_alphabeta_only_test.png") # store the output for later comparison bench = "mg_general_poisson_dirichlet" my_data = a.get_solution_object() if store_bench: my_data.write(f"{bench_dir}/{bench}") # do we do a comparison? if comp_bench: compare_file = f"{bench_dir}/{bench}" msg.warning("comparing to: %s " % (compare_file)) bench = io.read(compare_file) result = compare.compare(my_data, bench, rtol) if result == 0: msg.success(f"results match benchmark to within relative tolerance of {rtol}\n") else: msg.warning("ERROR: " + compare.errors[result] + "\n") return result # normal return -- error wrt true solution return enorm
[docs] def main(): N = [16, 32, 64] err = [] plot = False store = False do_compare = False for nx in N: if nx == max(N): plot = True enorm = test_general_poisson_dirichlet(nx, make_plot=plot, store_bench=store, comp_bench=do_compare) err.append(enorm) # plot the convergence N = np.array(N, dtype=np.float64) err = np.array(err) plt.clf() plt.loglog(N, err, "x", color="r") plt.loglog(N, err[0]*(N[0]/N)**2, "--", color="k") plt.xlabel("N") plt.ylabel("error") fig = plt.gcf() fig.set_size_inches(7.0, 6.0) plt.tight_layout() plt.savefig("mg_general_alphabeta_only_converge.png")
if __name__ == "__main__": main()