agd.ExportedCode.Notebooks_NonDiv.NonlinearMonotoneFirst2D
1# Code automatically exported from notebook NonlinearMonotoneFirst2D.ipynb in directory Notebooks_NonDiv 2# Do not modify 3from ... import Selling 4from ... import LinearParallel as lp 5from ... import AutomaticDifferentiation as ad 6from ... import Domain 7 8import numpy as np 9import matplotlib.pyplot as plt 10import scipy.linalg; 11import itertools 12 13def Gradient(u,A,bc,decomp=None): 14 """ 15 Approximates grad u(x), using finite differences along the axes of A. 16 """ 17 coefs,offsets = Selling.Decomposition(A) if decomp is None else decomp 18 du = bc.DiffCentered(u,offsets) 19 AGrad = lp.dot_AV(offsets.astype(float),(coefs*du)) # Approximates A * grad u 20 return lp.solve_AV(A,AGrad) # Approximates A^{-1} (A * grad u) = grad u 21 22def SchemeCentered(u,A,F,rhs,bc): 23 """ 24 Discretization of - Tr(A(x) hess u(x)) + F(grad u(x)) - rhs, 25 with Dirichlet boundary conditions. The scheme is second order, 26 and degenerate elliptic under suitable assumptions. 27 """ 28 # Compute the tensor decomposition 29 coefs,offsets = Selling.Decomposition(A) 30 A,coefs,offsets = (bc.as_field(e) for e in (A,coefs,offsets)) 31 32 # Obtain the first and second order finite differences 33 grad = Gradient(u,A,bc,decomp=(coefs,offsets)) 34 d2u = bc.Diff2(u,offsets) 35 36 # Numerical scheme in interior 37 residue = -lp.dot_VV(coefs,d2u) + F(grad) - rhs 38 39 # Placeholders outside domain 40 return np.where(bc.interior,residue,u-bc.grid_values) 41 42# Specialization for the quadratic non-linearity 43def SchemeCentered_Quad(u,A,omega,D,rhs,bc): 44 omega,D = (bc.as_field(e) for e in (omega,D)) 45 def F(g): return lp.dot_VAV(g-omega,D,g-omega) 46 return SchemeCentered(u,A,F,rhs,bc) 47 48def SchemeUpwind(u,A,omega,D,rhs,bc): 49 """ 50 Discretization of -Tr(A(x) hess u(x)) + \| grad u(x) - omega(x) \|_D(x)^2 - rhs, 51 with Dirichlet boundary conditions, using upwind finite differences for the first order part. 52 The scheme is degenerate elliptic if A and D are positive definite. 53 """ 54 # Compute the decompositions (here offset_e = offset_f) 55 nothing = (np.full((0,),0.), np.full((2,0),0)) # empty coefs and offsets 56 mu,offset_e = nothing if A is None else Selling.Decomposition(A) 57 nu,offset_f = nothing if D is None else Selling.Decomposition(D) 58 omega_f = lp.dot_VA(omega,offset_f.astype(float)) 59 60 # First and second order finite differences 61 maxi = np.maximum 62 mu,nu,omega_f = (bc.as_field(e) for e in (mu,nu,omega_f)) 63 64 dup = bc.DiffUpwind(u, offset_f) 65 dum = bc.DiffUpwind(u,-offset_f) 66 dup[...,bc.not_interior]=0. # Placeholder values to silence NaN warnings 67 dum[...,bc.not_interior]=0. 68 69 d2u = bc.Diff2(u,offset_e) 70 71 # Scheme in the interior 72 du = maxi(0.,maxi( omega_f - dup, -omega_f - dum) ) 73 residue = - lp.dot_VV(mu,d2u) + lp.dot_VV(nu,du**2) - rhs 74 75 # Placeholders outside domain 76 return np.where(bc.interior,residue,u-bc.grid_values)
def
Gradient(u, A, bc, decomp=None):
14def Gradient(u,A,bc,decomp=None): 15 """ 16 Approximates grad u(x), using finite differences along the axes of A. 17 """ 18 coefs,offsets = Selling.Decomposition(A) if decomp is None else decomp 19 du = bc.DiffCentered(u,offsets) 20 AGrad = lp.dot_AV(offsets.astype(float),(coefs*du)) # Approximates A * grad u 21 return lp.solve_AV(A,AGrad) # Approximates A^{-1} (A * grad u) = grad u
Approximates grad u(x), using finite differences along the axes of A.
def
SchemeCentered(u, A, F, rhs, bc):
23def SchemeCentered(u,A,F,rhs,bc): 24 """ 25 Discretization of - Tr(A(x) hess u(x)) + F(grad u(x)) - rhs, 26 with Dirichlet boundary conditions. The scheme is second order, 27 and degenerate elliptic under suitable assumptions. 28 """ 29 # Compute the tensor decomposition 30 coefs,offsets = Selling.Decomposition(A) 31 A,coefs,offsets = (bc.as_field(e) for e in (A,coefs,offsets)) 32 33 # Obtain the first and second order finite differences 34 grad = Gradient(u,A,bc,decomp=(coefs,offsets)) 35 d2u = bc.Diff2(u,offsets) 36 37 # Numerical scheme in interior 38 residue = -lp.dot_VV(coefs,d2u) + F(grad) - rhs 39 40 # Placeholders outside domain 41 return np.where(bc.interior,residue,u-bc.grid_values)
Discretization of - Tr(A(x) hess u(x)) + F(grad u(x)) - rhs, with Dirichlet boundary conditions. The scheme is second order, and degenerate elliptic under suitable assumptions.
def
SchemeCentered_Quad(u, A, omega, D, rhs, bc):
def
SchemeUpwind(u, A, omega, D, rhs, bc):
49def SchemeUpwind(u,A,omega,D,rhs,bc): 50 """ 51 Discretization of -Tr(A(x) hess u(x)) + \| grad u(x) - omega(x) \|_D(x)^2 - rhs, 52 with Dirichlet boundary conditions, using upwind finite differences for the first order part. 53 The scheme is degenerate elliptic if A and D are positive definite. 54 """ 55 # Compute the decompositions (here offset_e = offset_f) 56 nothing = (np.full((0,),0.), np.full((2,0),0)) # empty coefs and offsets 57 mu,offset_e = nothing if A is None else Selling.Decomposition(A) 58 nu,offset_f = nothing if D is None else Selling.Decomposition(D) 59 omega_f = lp.dot_VA(omega,offset_f.astype(float)) 60 61 # First and second order finite differences 62 maxi = np.maximum 63 mu,nu,omega_f = (bc.as_field(e) for e in (mu,nu,omega_f)) 64 65 dup = bc.DiffUpwind(u, offset_f) 66 dum = bc.DiffUpwind(u,-offset_f) 67 dup[...,bc.not_interior]=0. # Placeholder values to silence NaN warnings 68 dum[...,bc.not_interior]=0. 69 70 d2u = bc.Diff2(u,offset_e) 71 72 # Scheme in the interior 73 du = maxi(0.,maxi( omega_f - dup, -omega_f - dum) ) 74 residue = - lp.dot_VV(mu,d2u) + lp.dot_VV(nu,du**2) - rhs 75 76 # Placeholders outside domain 77 return np.where(bc.interior,residue,u-bc.grid_values)
Discretization of -Tr(A(x) hess u(x)) + \| grad u(x) - omega(x) \|_D(x)^2 - rhs, with Dirichlet boundary conditions, using upwind finite differences for the first order part. The scheme is degenerate elliptic if A and D are positive definite.