<jupyter_start><jupyter_text>3. 1D Burgers EquationWe consider the 1d Burgers equation$$\partial_t u + u \partial_x u = \nu \frac{\partial ^2u}{\partial x^2}$$$u_0(x) := u(x,t)$ denotes the initial condition.We choose homogeneous neumann boundary conditions in this example, i.e.$$\partial_n u = 0, \partial \Omega$$ with $\Omega = (0,1)$<jupyter_code># needed imports
from numpy import zeros, ones, linspace, zeros_like
from matplotlib.pyplot import plot, show
%matplotlib inline
# Initial condition
from numpy import exp
#u0 = lambda x: exp(-(x-.3)**2/.05**2)
u0 = lambda x: exp(-(x-.5)**2/.02**2)

grid = linspace(0., 1., 401)
u = u0(grid)
plot(grid, u) ; show()<jupyter_output><empty_output><jupyter_text>Time schemeWe shall use a $\theta$-scheme in this case and consider the following problem$$\frac{u^{n+1}-u^n}{\Delta t} + \theta~ u^{n+1} \partial_x u^{n+1} + (1-\theta)~ u^n \partial_x u^n = \theta~\nu \frac{\partial ^2u^{n+1}}{\partial x^2} + (1-\theta)~\nu \frac{\partial ^2u^{n}}{\partial x^2}$$hence$$u^{n+1} + \Delta t ~ \theta~ u^{n+1} \partial_x u^{n+1} - \Delta t ~ \theta~\nu \frac{\partial ^2u^{n+1}}{\partial x^2} =u^{n} - \Delta t ~ (1-\theta)~ u^{n} \partial_x u^{n} + \Delta t ~ (1-\theta)~\nu \frac{\partial ^2u^{n}}{\partial x^2}$$from now on, we shall denote by $f^n$ the right hand side of the previous equation$$f^n := u^{n} - \Delta t ~ (1-\theta)~ u^{n} \partial_x u^{n} + \Delta t ~ (1-\theta)~\nu \frac{\partial ^2u^{n}}{\partial x^2}$$ Weak formulationLet $v \in \mathcal{V}$ be a function test, we have by integrating by parts the highest order term:$$\langle v, u^{n+1}  \rangle + \Delta t ~ \theta~ \langle v, u^{n+1} \partial_x u^{n+1}  \rangle+ \Delta t ~ \theta~\nu \langle \frac{\partial v}{\partial x}, \frac{\partial u^{n+1}}{\partial x}  \rangle=\langle v, f^n \rangle$$The previous weak formulation is still nonlinear with respect to $u^{n+1}$. We shall then follow the same strategy as for the previous chapter on nonlinear Poisson problem.The strategy is to define the left hand side as a **LinearForm** with respect to $v$, then linearize it around $u^{n+1}$. We therefor can use either Picard or Newton method to treat the nonlinearity.We consider the following linear form$$G(v;u,w) := \langle v, u  \rangle + \Delta t ~ \theta~ \langle v, w \partial_x u  \rangle+ \Delta t ~ \theta~\nu \langle \frac{\partial v}{\partial x}, \frac{\partial u}{\partial x}  \rangle, \quad \forall u,v,w \in \mathcal{V}$$Our problem is then$$\mbox{Find } u^{n+1} \in \mathcal{V}, \mbox{such that}\\G(v;u^{n+1},u^{n+1}) = l(v), \quad \forall v \in \mathcal{V}$$where$$l(v) := \int_{\Omega} f^n v ~d\Omega, \quad \forall v \in \mathcal{V}$$  Abstract Model<jupyter_code>from sympde.core import Constant
from sympde.expr import BilinearForm, LinearForm, integral
from sympde.topology import ScalarFunctionSpace, Line, element_of
from sympde.topology import dx1 as dx # TODO: this is a bug right now
from sympde.expr     import find
from sympde.expr.expr import linearize

from psydac.fem.basic          import FemField
from psydac.api.discretization import discretize
domain = Line()

V = ScalarFunctionSpace('V', domain)
u  = element_of(V, name='u')
v  = element_of(V, name='v')
w  = element_of(V, name='w')
un  = element_of(V, name='un') # time iteration
uk  = element_of(V, name='uk') # nonlinear solver iteration

x = domain.coordinates

nu = Constant('nu')
theta = Constant('theta')
dt = Constant('dt')<jupyter_output><empty_output><jupyter_text>Defining the Linear form $G$<jupyter_code># Linear form g: V --> R
expr = v * u + dt*theta*v*w*dx(u) + dt*theta*nu*dx(v)*dx(u)
g = LinearForm(v, integral(domain, expr))<jupyter_output><empty_output><jupyter_text>Defining the Linear form $l$<jupyter_code># Linear form l: V --> R
expr = v * un - dt*theta*v*un*dx(un) - dt*theta*nu*dx(v)*dx(un)
l = LinearForm(v, integral(domain, expr))<jupyter_output><empty_output><jupyter_text>Picard Method$$\mbox{Find } u_{n+1} \in \mathcal{V}_h, \mbox{such that}\\G(v;u_{n+1},u_n) = l(v), \quad \forall v \in \mathcal{V}_h$$  Picard iteration<jupyter_code># Variational problem
picard = find(u, forall=v, lhs=g(v, u=u,w=uk), rhs=l(v))<jupyter_output><empty_output><jupyter_text>Newton MethodLet's define $$F(v;u) := G(v;u,u) -l(v), \quad \forall v \in \mathcal{V}$$Newton method writes$$\mbox{Find } u_{n+1} \in \mathcal{V}_h, \mbox{such that}\\F^{\prime}(\delta u,v; u_n) = - F(v;u_n), \quad \forall v \in \mathcal{V} \\u_{n+1} := u_{n} + \delta u, \quad \delta u \in \mathcal{V}$$  Newton Iteration<jupyter_code>F = LinearForm(v, g(v,w=u)-l(v))
du  = element_of(V, name='du')

Fprime = linearize(F, u, trials=du)

# Variational problem
newton = find(du, forall=v, lhs=Fprime(du, v,u=uk), rhs=-F(v,u=uk))<jupyter_output><empty_output><jupyter_text>Discrete Space<jupyter_code># Create computational domain from topological domain
domain_h = discretize(domain, ncells=[64], comm=None)

# Discrete spaces
Vh = discretize(V, domain_h, degree=[2])<jupyter_output><empty_output><jupyter_text>Picard method<jupyter_code># Discretize equation 
picard_h = discretize(picard, domain_h, [Vh, Vh])

def picard(Un, niter=10):
    Uk = FemField( Vh, Vh.vector_space.zeros() )

    for i in range(niter):
        Uk  = picard_h.solve(uk=Uk, un=Un)
        
    return Uk<jupyter_output><empty_output><jupyter_text>Newton method<jupyter_code># Discretize equation 
newton_h = discretize(newton, domain_h, [Vh, Vh])

def newton(Un, niter=10):
    Uk = FemField( Vh, Vh.vector_space.zeros() )

    for i in range(niter):
        delta_x  = newton_h.solve(uk=Uk, un=Un)
        Uk = FemField( Vh, delta_x.coeffs + Uk.coeffs )
        
    return Uk<jupyter_output><empty_output><jupyter_text>L2-Projection<jupyter_code>a_mass = BilinearForm((u,v), integral(domain , u*v))

from sympy import exp
#u0 = exp(-(x-.3)**2/.05**2)
u0 = exp(-(x-.5)**2/.02**2)

l_mass = LinearForm(v, integral(domain, u0*v))

# Abstract projection
projection = find(u, forall=v, lhs=a_mass(u, v), rhs=l_mass(v))

# Discrete projection
projection_h = discretize(projection, domain_h, [Vh, Vh])
u0_h = projection_h.solve()
from utilities.plot import plot_field_1d
plot_field_1d(Vh.knots[0], Vh.degree[0], u0_h.coeffs.toarray(), nx=401)<jupyter_output><empty_output><jupyter_text>Important noteRight now, Psydac does not provide a **grmes** solver. For this reason, we shall rewrite our nonlinear solver using the following **gmres_driver**<jupyter_code>from scipy.sparse.linalg import gmres
def gmres_driver(equation_h, tol=1.e-8, maxiter=5000):
    M   = equation_h.linear_system.lhs.tosparse()
    rhs = equation_h.linear_system.rhs.toarray()
    
    x, status = gmres(M, rhs, tol=tol, maxiter=maxiter)
    
    xh = FemField( Vh, Vh.vector_space.zeros() )
    n = xh.coeffs._data.shape[0]
    xh.coeffs._data[Vh.degree[0]:n-Vh.degree[0]] = x[:]
    
    return xh
def picard(Un, theta, nu, dt, niter=10):
    Uk = FemField( Vh, Vh.vector_space.zeros() )

    for i in range(niter):
        picard_h.assemble(uk=Uk, un=Un, theta=theta, nu=nu, dt=dt)
        Uk = gmres_driver(picard_h, tol=1.e-8, maxiter=5000)
        
    return Uk
def newton(Un, theta, nu, dt, niter=10):
    Uk = FemField( Vh, Vh.vector_space.zeros() )

    for i in range(niter):
        newton_h.assemble(uk=Uk, un=Un, theta=theta, nu=nu, dt=dt)
        delta_x = gmres_driver(newton_h, tol=1.e-8, maxiter=5000)
        Uk = FemField( Vh, delta_x.coeffs + Uk.coeffs )
        
    return Uk
Uk = picard(u0_h, theta=0.5, nu=0.5, dt=0.01)
plot_field_1d(Vh.knots[0], Vh.degree[0], Uk.coeffs.toarray(), nx=401)
Uk = newton(u0_h, theta=0.5, nu=0.5, dt=0.01)
plot_field_1d(Vh.knots[0], Vh.degree[0], Uk.coeffs.toarray(), nx=401)
def time_solver(theta, nu, dt, T, nonlinear_solver, U0):
    n_time = int(T / dt)
    Uk = U0.copy()
    for i_time in range(0, n_time):
        Uk = nonlinear_solver(Uk, theta=theta, nu=nu, dt=dt)
    return Uk
Uk = time_solver(theta=0.5, nu=0.02, dt=0.01, T=0.1, nonlinear_solver=newton, U0=u0_h)
plot_field_1d(Vh.knots[0], Vh.degree[0], Uk.coeffs.toarray(), nx=401)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1734820322180917, 'fraction_numerical': 0.02094175960346964, 'mean_word_length': 4.550894085281981, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '528593', 'license': 'MIT', 'n_tokens_mistral': 3341, 'n_tokens_neox': 3198, 'n_words': 778, 'path': 'lessons/Chapter3/02_burgers_1d.ipynb', 'repo_name': 'pyccel/IGA-Python', 'size': '8070'}