<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Multigrid examples
Step1:  Constant-coefficent Poisson equation
We want to solve 
$$\phi_{xx} + \phi_{yy} = -2[(1-6x^2)y^2(1-y^2) + (1-6y^2)x^2(1-x^2)]$$
on 
$$[0,1]\times [0,1]$$ 
with homogeneous Dirichlet boundary conditions (this example comes from "A Multigrid Tutorial").
This has the analytic solution
$$u(x,y) = (x^2 - x^4)(y^4 - y^2)$$
We start by setting up a multigrid object--this needs to know the number of zones our problem is defined on
Step2:  Next, we initialize the RHS.  To make life easier, the CellCenterMG2d object has the coordinates of the solution grid (including ghost cells) as mg.x2d and mg.y2d (these are two-dimensional arrays).
Step3:  The last setup step is to initialize the solution--this is the starting point for the solve.  Usually we just want to start with all zeros, so we use the init_zeros() method
Step4:  we can now solve -- there are actually two different techniques we can do here.  We can just do pure smoothing on the solution grid using mg.smooth(mg.nlevels-1, N), where N is the number of smoothing iterations.  To get the solution N will need to be large and this will take a long time.
Multigrid accelerates the smoothing.  We can do a V-cycle multigrid solution using mg.solve()
Step5:  We can access the solution on the finest grid using get_solution()
Step6:  we can also get the gradient of the solution
Step7:  General linear elliptic equation
The GeneralMG2d class implements support for a general elliptic equation of the form
Step8:  For reference, we'll define a function providing the analytic solution
Step9:  Now the coefficents--note that since $\gamma$ is a vector, we have a different function for each component
Step10:  and the righthand side function
Step11:  Our inhomogeneous boundary conditions require a function that can be evaluated on the boundary to give the value
Step12:  Now we can setup our grid object and the coefficients, which are stored as a CellCenter2d object.  Note, the coefficients do not need to have the same boundary conditions as $\phi$ (and for real problems, they may not).  The one that matters the most is $\beta$, since that will need to be averaged to the edges of the domain, so the boundary conditions on the coefficients are important.
Here we use Neumann boundary conditions
Step13:  Now we can setup the multigrid object
Step14:  just as before, we specify the righthand side and initialize the solution
Step15:  and we can solve it
Step16:  We can compare to the true solution
Step17:  The norm of the error is
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
from __future__ import print_function
import numpy as np
import mesh.boundary as bnd
import mesh.patch as patch
import multigrid.MG as MG
Explanation: Multigrid examples
End of explanation
nx = ny = 256
mg = MG.CellCenterMG2d(nx, ny, 
                       xl_BC_type="dirichlet", xr_BC_type="dirichlet",
                       yl_BC_type="dirichlet", yr_BC_type="dirichlet", verbose=1)
Explanation: Constant-coefficent Poisson equation
We want to solve 
$$\phi_{xx} + \phi_{yy} = -2[(1-6x^2)y^2(1-y^2) + (1-6y^2)x^2(1-x^2)]$$
on 
$$[0,1]\times [0,1]$$ 
with homogeneous Dirichlet boundary conditions (this example comes from "A Multigrid Tutorial").
This has the analytic solution
$$u(x,y) = (x^2 - x^4)(y^4 - y^2)$$
We start by setting up a multigrid object--this needs to know the number of zones our problem is defined on
End of explanation
def rhs(x, y):
    return -2.0*((1.0-6.0*x**2)*y**2*(1.0-y**2) + (1.0-6.0*y**2)*x**2*(1.0-x**2))
mg.init_RHS(rhs(mg.x2d, mg.y2d))
Explanation: Next, we initialize the RHS.  To make life easier, the CellCenterMG2d object has the coordinates of the solution grid (including ghost cells) as mg.x2d and mg.y2d (these are two-dimensional arrays).
End of explanation
mg.init_zeros()
Explanation: The last setup step is to initialize the solution--this is the starting point for the solve.  Usually we just want to start with all zeros, so we use the init_zeros() method
End of explanation
mg.solve()
Explanation: we can now solve -- there are actually two different techniques we can do here.  We can just do pure smoothing on the solution grid using mg.smooth(mg.nlevels-1, N), where N is the number of smoothing iterations.  To get the solution N will need to be large and this will take a long time.
Multigrid accelerates the smoothing.  We can do a V-cycle multigrid solution using mg.solve()
End of explanation
phi = mg.get_solution()
plt.imshow(np.transpose(phi.v()), origin="lower")
Explanation: We can access the solution on the finest grid using get_solution()
End of explanation
gx, gy = mg.get_solution_gradient()
plt.subplot(121)
plt.imshow(np.transpose(gx.v()), origin="lower")
plt.subplot(122)
plt.imshow(np.transpose(gy.v()), origin="lower")
Explanation: we can also get the gradient of the solution
End of explanation
import multigrid.general_MG as gMG
Explanation: General linear elliptic equation
The GeneralMG2d class implements support for a general elliptic equation of the form:
$$\alpha \phi + \nabla \cdot (\beta \nabla \phi) + \gamma \cdot \nabla \phi = f$$
with inhomogeneous boundary condtions.
It subclasses the CellCenterMG2d class, and the basic interface is the same
We will solve the above with
\begin{align}
\alpha &= 10 \
\beta &= xy + 1 \
\gamma &= \hat{x} + \hat{y}
\end{align}
and
\begin{align}
f = &-\frac{\pi}{2}(x + 1)\sin\left(\frac{\pi y}{2}\right) \cos\left(\frac{\pi x}{2}\right ) \
    &-\frac{\pi}{2}(y + 1)\sin\left(\frac{\pi x}{2}\right) \cos\left(\frac{\pi y}{2}\right ) \
    &+\left(\frac{-\pi^2 (xy+1)}{2} + 10\right) \cos\left(\frac{\pi x}{2}\right) \cos\left(\frac{\pi y}{2}\right)
\end{align}
on $[0, 1] \times [0,1]$ with boundary conditions:
\begin{align}
\phi(x=0) &= \cos(\pi y/2) \
\phi(x=1) &= 0 \
\phi(y=0) &= \cos(\pi x/2) \
\phi(y=1) &= 0
\end{align}
This has the exact solution:
$$\phi = \cos(\pi x/2) \cos(\pi y/2)$$
End of explanation
def true(x,y):                                                                                         
    return np.cos(np.pi*x/2.0)*np.cos(np.pi*y/2.0)                                                     
Explanation: For reference, we'll define a function providing the analytic solution
End of explanation
def alpha(x,y):                                                                                        
    return 10.0*np.ones_like(x)                                                                        
                                                                                                       
def beta(x,y):                                                                                         
    return x*y + 1.0                                                                                   
                                                                                                       
def gamma_x(x,y):                                                                                      
    return np.ones_like(x)                                                                             
                                                                                                       
def gamma_y(x,y):                                                                                      
    return np.ones_like(x)                                                                             
Explanation: Now the coefficents--note that since $\gamma$ is a vector, we have a different function for each component
End of explanation
def f(x,y):
    return -0.5*np.pi*(x + 1.0)*np.sin(np.pi*y/2.0)*np.cos(np.pi*x/2.0) - \
            0.5*np.pi*(y + 1.0)*np.sin(np.pi*x/2.0)*np.cos(np.pi*y/2.0) + \
            (-np.pi**2*(x*y+1.0)/2.0 + 10.0)*np.cos(np.pi*x/2.0)*np.cos(np.pi*y/2.0)
Explanation: and the righthand side function
End of explanation
def xl_func(y):
    return np.cos(np.pi*y/2.0)
def yl_func(x):
    return np.cos(np.pi*x/2.0)
Explanation: Our inhomogeneous boundary conditions require a function that can be evaluated on the boundary to give the value
End of explanation
import mesh.patch as patch
nx = ny = 128
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)
Explanation: Now we can setup our grid object and the coefficients, which are stored as a CellCenter2d object.  Note, the coefficients do not need to have the same boundary conditions as $\phi$ (and for real problems, they may not).  The one that matters the most is $\beta$, since that will need to be averaged to the edges of the domain, so the boundary conditions on the coefficients are important.
Here we use Neumann boundary conditions
End of explanation
a = gMG.GeneralMG2d(nx, ny,
                    xl_BC_type="dirichlet", yl_BC_type="dirichlet",
                    xr_BC_type="dirichlet", yr_BC_type="dirichlet",
                    xl_BC=xl_func,
                    yl_BC=yl_func,
                    coeffs=d,
                    verbose=1, vis=0, true_function=true)
Explanation: Now we can setup the multigrid object
End of explanation
a.init_zeros()
a.init_RHS(f(a.x2d, a.y2d))
Explanation: just as before, we specify the righthand side and initialize the solution
End of explanation
a.solve(rtol=1.e-10)
Explanation: and we can solve it
End of explanation
v = a.get_solution()
b = true(a.x2d, a.y2d)
e = v - b
Explanation: We can compare to the true solution
End of explanation
e.norm()
Explanation: The norm of the error is
End of explanation
<END_TASK>