<jupyter_start><jupyter_text>Non-Linear Poisson in 2D In this section, we consider the non-linear Poisson problem:$$-\nabla \cdot \left( (1+u^2) \nabla u \right) = f, \Omega\\u = 0, \partial \Omega$$where $\Omega$ denotes the unit square.For testing, we shall take a function $u$ that fulfills the boundary condition, the compute $f$ as$$f(x,y) = -\nabla^2 u - F(u)$$ The weak formulation is$$\int_{\Omega} (1+u^2) \nabla u \cdot \nabla v ~ d\Omega = \int_{\Omega} f v ~d\Omega, \quad \forall v \in \mathcal{V}$$ For the sack of generality, we shall consider the linear form$$G(v;u,w) := \int_{\Omega} (1+w^2) \nabla u \cdot \nabla v ~ d\Omega, \quad \forall u,v,w \in \mathcal{V}$$ Our problem is then$$\mbox{Find } u \in \mathcal{V}, \mbox{such that}\\G(v;u,u) = l(v), \quad \forall v \in \mathcal{V}$$where$$l(v) := \int_{\Omega} f v ~d\Omega, \quad \forall v \in \mathcal{V}$$  Topological domain<jupyter_code>domain = Square()
B_dirichlet_0 = domain.boundary<jupyter_output><empty_output><jupyter_text>Function Space<jupyter_code>V  = ScalarFunctionSpace('V', domain)<jupyter_output><empty_output><jupyter_text>Defining the Linear form $G$<jupyter_code>u  = element_of(V, name='u')
v  = element_of(V, name='v')
w  = element_of(V, name='w')

# Linear form g: V --> R
g = LinearForm(v, integral(domain, (1+w**2)*dot(grad(u), grad(v))))<jupyter_output><empty_output><jupyter_text>Defining the Linear form L<jupyter_code>solution = sin(pi*x)*sin(pi*y)
f = 2*pi**2*(sin(pi*x)**2*sin(pi*y)**2 + 1)*sin(pi*x)*sin(pi*y) - 2*pi**2*sin(pi*x)**3*sin(pi*y)*cos(pi*y)**2 - 2*pi**2*sin(pi*x)*sin(pi*y)**3*cos(pi*x)**2
# Linear form l: V --> R
l = LinearForm(v, integral(domain, f * v))<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$$  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}$$  Computing $F^{\prime}$ the derivative of $F$**SymPDE** allows you to linearize a linear form and get a bilinear form, using the function **linearize**<jupyter_code>F = LinearForm(v, g(v,w=u)-l(v))
du  = element_of(V, name='du')

Fprime = linearize(F, u, trials=du)<jupyter_output><empty_output><jupyter_text>Picard Method  Abstract Model<jupyter_code>un  = element_of(V, name='un')

# Bilinear form a: V x V --> R
a = BilinearForm((u, v), g(v, u=u,w=un))

# Dirichlet boundary conditions
bc = [EssentialBC(u, 0, B_dirichlet_0)]

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)

# Error norms
error  = u - solution
l2norm = Norm(error, domain, kind='l2')<jupyter_output><empty_output><jupyter_text>Discretization<jupyter_code># Create computational domain from topological domain
domain_h = discretize(domain, ncells=[16,16], comm=None)

# Discrete spaces
Vh = discretize(V, domain_h, degree=[2,2])

# Discretize equation using Dirichlet bc
equation_h = discretize(equation, domain_h, [Vh, Vh])

# Discretize error norms
l2norm_h = discretize(l2norm, domain_h, Vh)<jupyter_output><empty_output><jupyter_text>Picard solver<jupyter_code>def picard(niter=10):
    Un = FemField( Vh, Vh.vector_space.zeros() )

    for i in range(niter):
        Un  = equation_h.solve(un=Un)

        # Compute error norms
        l2_error = l2norm_h.assemble(u=Un)

        print('l2_error = ', l2_error)
        
    return Un
Un = picard(niter=5)
from matplotlib import pyplot as plt
from utilities.plot import plot_field_2d

nbasis = [w.nbasis for w in Vh.spaces]
p1,p2 = Vh.degree
x = Un.coeffs._data[p1:-p1,p2:-p2]
u = x.reshape(nbasis)
plot_field_2d(Vh.knots, Vh.degree, u) ; plt.colorbar()<jupyter_output><empty_output><jupyter_text>Newton Method  Abstract Model<jupyter_code># Dirichlet boundary conditions
bc = [EssentialBC(du, 0, B_dirichlet_0)]

# Variational problem
equation   = find(du, forall=v, lhs=Fprime(du, v,u=un), rhs=-F(v,u=un), bc=bc)<jupyter_output><empty_output><jupyter_text>Discretization<jupyter_code># Create computational domain from topological domain
domain_h = discretize(domain, ncells=[16,16], comm=None)

# Discrete spaces
Vh = discretize(V, domain_h, degree=[2,2])

# Discretize equation using Dirichlet bc
equation_h = discretize(equation, domain_h, [Vh, Vh])

# Discretize error norms
l2norm_h = discretize(l2norm, domain_h, Vh)<jupyter_output><empty_output><jupyter_text>Newton Solver<jupyter_code>def newton(niter=10):
    Un = FemField( Vh, Vh.vector_space.zeros() )

    for i in range(niter):
        delta_x  = equation_h.solve(un=Un)
        Un = FemField( Vh, delta_x.coeffs + Un.coeffs )

        # Compute error norms
        l2_error = l2norm_h.assemble(u=Un)

        print('l2_error = ', l2_error)
        
    return Un
un = newton(niter=5)
nbasis = [w.nbasis for w in Vh.spaces]
p1,p2 = Vh.degree
x = un.coeffs._data[p1:-p1,p2:-p2]
u = x.reshape(nbasis)
plot_field_2d(Vh.knots, Vh.degree, u) ; plt.colorbar()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16911764705882354, 'fraction_numerical': 0.015479876160990712, 'mean_word_length': 5.05978898007034, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, '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': '9368', 'license': 'MIT', 'n_tokens_mistral': 2161, 'n_tokens_neox': 2066, 'n_words': 464, 'path': 'lessons/Chapter3/01_nonlinear_poisson_2d.ipynb', 'repo_name': 'pyccel/IGA-Python', 'size': '5168'}