agd.ExportedCode.Notebooks_NonDiv.LinearMonotoneSchemes2D

18def norm(arr,ord=2,axis=None,keepdims=False,averaged=False):
19	"""
20	Returns L^p norm of array, seen as a vector, w.r.t. weights.
21	Defined as : (sum_i x[i]^p)^(1/p)
22
23	Remark : not a matrix operator norm
24	
25	Inputs:
26	 - ord : exponent p
27	 - axis : int or None, axis along which to compute the norm. 
28	 - keepdims : wether to keep singleton dimensions.
29	 - averaged : wether to introduce a normalization factor, so that norm(ones(...))=1
30
31	Compatible with automatic differentiation.
32	"""
33	arr = ad_generic.array(arr)
34	if ord==np.inf or ord%2!=0:
35		try: arr = np.abs(arr)
36		except TypeError: arr = np.vectorize(np.abs)(arr)
37	if ord==np.inf: return np.max(arr,axis=axis,keepdims=keepdims)
38	sum_pow = np.sum(arr**ord,axis=axis,keepdims=keepdims)
39	
40	if averaged:
41		size = arr.size if axis is None else arr.shape[axis]
42		sum_pow/=size
43
44	return sum_pow**(1./ord)

15def streamplot_ij(X,Y,VX,VY,subsampling=1,*varargs,**kwargs):
16    def f(array): return array[::subsampling,::subsampling].T
17    return plt.streamplot(f(X),f(Y),f(VX),f(VY),*varargs,**kwargs) # Transpose everything

19def SchemeCentered(u,cst,mult,omega,diff,bc,ret_hmax=False):
20    """Discretization of a linear non-divergence form second order PDE
21        cst + mult u + <omega,grad u>- tr(diff hess(u)) = 0
22        Second order accurate, centered yet monotone finite differences are used for <omega,grad u>
23        - bc : boundary conditions. 
24        - ret_hmax : return the largest grid scale for which monotony holds
25    """
26    # Decompose the tensor field
27    coefs2,offsets = Selling.Decomposition(diff)
28    
29    # Decompose the vector field
30    scals = lp.dot_VA(lp.solve_AV(diff,omega), offsets.astype(float))
31    coefs1 = coefs2*scals
32    if ret_hmax: return 2./norm(scals,ord=np.inf)
33    
34    # Compute the first and second order finite differences    
35    du  = bc.DiffCentered(u,offsets)
36    d2u = bc.Diff2(u,offsets)
37    
38    # In interior : cst + mult u + <omega,grad u>- tr(diff hess(u)) = 0
39    coefs1,coefs2 = (bc.as_field(e) for e in (coefs1,coefs2))    
40    residue = cst + mult*u +lp.dot_VV(coefs1,du) - lp.dot_VV(coefs2,d2u)
41    
42    # On boundary : u-bc = 0
43    return np.where(bc.interior,residue,u-bc.grid_values)

45def SchemeUpwind(u,cst,mult,omega,diff,bc):
46    """Discretization of a linear non-divergence form second order PDE
47        cst + mult u + <omega,grad u>- tr(diff hess(u)) = 0
48        First order accurate, upwind finite differences are used for <omega,grad u>
49        - bc : boundary conditions. 
50    """
51    # Decompose the tensor field
52    coefs2,offsets2 = Selling.Decomposition(diff)
53    omega,coefs2 = (bc.as_field(e) for e in (omega,coefs2))    
54
55    # Decompose the vector field
56    coefs1 = -np.abs(omega)
57    basis = bc.as_field(np.eye(len(omega)))
58    offsets1 = -np.sign(omega)*basis
59    
60    # Compute the first and second order finite differences    
61    du  = bc.DiffUpwind(u,offsets1.astype(int))
62    d2u = bc.Diff2(u,offsets2)
63    
64    # In interior : cst + mult u + <omega,grad u>- tr(diff hess(u)) = 0
65    residue = cst + mult*u +lp.dot_VV(coefs1,du) - lp.dot_VV(coefs2,d2u)
66    
67    # On boundary : u-bc = 0
68    return np.where(bc.interior,residue,u-bc.grid_values)