agd.LinearPDE

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4"""
  5This module is DEPRECATED, but kept for compatibility.
  6"""
  7
  8import numpy as np
  9import scipy.sparse as sp
 10
 11from . import Selling
 12from . import LinearParallel as LP
 13from .AutomaticDifferentiation.cupy_generic import get_array_module
 14
 15def OperatorMatrix(diff,omega=None,mult=None, \
 16		gridScale=1,
 17		boundaryConditions='Periodic', 
 18		divergenceForm=False,
 19		intrinsicDrift=False):
 20	r"""
 21	Constructs a linear operator sparse matrix, given as input 
 22	- an array `diff` of symmetric positive definite matrices, 
 23		with shape $(d,d,n_1,...,n_d)$ where $d$ is the domain dimension.
 24	- an array `omega` of vectors (optionnal), with shape $(d,n_1,...,n_d)$.
 25	- an array of scalars (optionnal), with shape $(n_1,...,n_d)$.
 26
 27	additional parameters
 28	- a grid scale 
 29	- boundary conditions, possibly axis by axis 
 30		('Periodic', 'Reflected', 'Neumann', 'Dirichlet') 
 31	- divergence form or not
 32
 33	The discretized operator is
 34	$$
 35	 {-} \mathrm{div}(D \nabla u) + < \\omega, \nabla u> + mult*u,
 36	$$
 37	denoting $D:=$`diff` and $\\omega:=$`omega`.
 38
 39	Replace the first term with $\mathrm{Tr}(D \nabla^2 u)$ in the 
 40	non-divergence form case.
 41	
 42	Returns : a list of triplets, for building a coo matrix
 43	"""
 44	# ----- Get the domain shape -----
 45	bounds = diff.shape[2:]
 46	dim = len(bounds)
 47	xp = get_array_module(diff)
 48	if isinstance(boundaryConditions,str):
 49		boundaryConditions = np.full( (dim,2), boundaryConditions)
 50	elif len(boundaryConditions)!=dim:
 51		raise ValueError("""OperatorMatrix error : 
 52		inconsistent boundary conditions""")
 53	
 54	if diff.shape[:2]!=(dim,dim):
 55		raise ValueError("OperatorMatrix error : inconsistent matrix dimensions")
 56
 57	# -------- Decompose the tensors --------
 58	coef,offset = Selling.Decomposition(diff)
 59	nCoef = coef.shape[0]
 60	
 61	# ------ Check bounds or apply periodic boundary conditions -------
 62	grid = np.mgrid[tuple(slice(0,n) for n in bounds)]
 63	grid = grid[:,None,...]
 64	grid = xp.asarray(grid,dtype=offset.dtype)
 65
 66	neighPos = grid + offset
 67	neighNeg = grid - offset
 68	neumannPos = np.full_like(coef,False,dtype=bool)
 69	neumannNeg = np.full_like(coef,False,dtype=bool)
 70	dirichletPos = np.full_like(coef,False,dtype=bool)
 71	dirichletNeg = np.full_like(coef,False,dtype=bool)
 72	
 73	for neigh_,neumann,dirichlet in zip( (neighPos,neighNeg), (neumannPos,neumannNeg), (dirichletPos,dirichletNeg) ): 
 74		for neigh,cond_,bound in zip(neigh_,boundaryConditions,bounds): # Component by component
 75			for out,cond in zip( (neigh<0,neigh>=bound), cond_):
 76				if cond=='Periodic':
 77   				 	neigh[out] %= bound 
 78				elif cond=='Neumann':
 79					neumann[out] = True
 80				elif cond=='Dirichlet':
 81					dirichlet[out] = True
 82	
 83	# ------- Get the neighbor indices --------
 84	# Cumulative product in reverse order, omitting last term, beginning with 1
 85	cum = tuple(list(reversed(np.cumprod(list(reversed(bounds+(1,)))))))[1:]
 86	bCum = np.broadcast_to( np.reshape(cum, (dim,)+(1,)*(dim+1)), offset.shape)
 87	bCum = xp.asarray(bCum,dtype=grid.dtype)
 88
 89	index = (grid*bCum).sum(0)
 90	indexPos = (neighPos*bCum).sum(0)
 91	indexNeg = (neighNeg*bCum).sum(0)
 92	index = np.broadcast_to(index,indexPos.shape)
 93	
 94	# ------- Get the coefficients for the first order term -----
 95	if omega is not None:
 96		if intrinsicDrift:
 97			eta=omega
 98		else:
 99			eta = LP.dot_AV(LP.inverse(diff),omega)
100			
101		scalEta = LP.dot_VV(offset.astype(float), 
102			np.broadcast_to(np.reshape(eta,(dim,1,)+bounds),offset.shape)) 
103		coefOmega = coef*scalEta
104
105	# ------- Create the triplets ------
106	
107	# Second order part
108	# Nemann : remove all differences which are not inside (a.k.a multiply coef by inside)
109	# TODO : Dirichlet : set to zero the coef only for the outside part
110	
111	coef = coef.reshape(-1)/ (gridScale**2) # Take grid scale into account
112
113	index = index.reshape(-1)
114	indexPos = indexPos.reshape(-1)
115	indexNeg = indexNeg.reshape(-1)
116
117	nff = lambda t : np.logical_not(t).astype(float).reshape(-1)
118	IP = nff(np.logical_or(neumannPos,dirichletPos))
119	IN = nff(np.logical_or(neumannNeg,dirichletNeg))
120	iP = nff(neumannPos)
121	iN = nff(neumannNeg)
122
123	if divergenceForm:
124		row = np.concatenate((index, indexPos, index, indexPos))
125		col = np.concatenate((index, index, indexPos, indexPos))
126		data = np.concatenate((iP*coef/2, -IP*coef/2, -IP*coef/2, IP*coef/2))
127		
128		row  = np.concatenate(( row, index, indexNeg, index, indexNeg))
129		col  = np.concatenate(( col, index, index, indexNeg, indexNeg))
130		data = np.concatenate((data, iN*coef/2, -IN*coef/2, -IN*coef/2, IN*coef/2))
131		
132	else:
133		row = np.concatenate( (index, index,	index))
134		col = np.concatenate( (index, indexPos, indexNeg))
135		data = np.concatenate((iP*coef+iN*coef, -IP*coef, -IN*coef))
136	
137
138	# First order part, using centered finite differences
139	if omega is not None:	   
140		coefOmega = coefOmega.flatten() / gridScale # Take grid scale in
141		row = np.concatenate((row, index,	index))
142		col = np.concatenate((col, indexPos, indexNeg))
143		data= np.concatenate((data,IP*iN*coefOmega/2,-IN*iP*coefOmega/2))
144	
145	if mult is not None:
146		# TODO Non periodic boundary conditions
147		size=np.prod(bounds)
148		row = np.concatenate((row, range(size)))
149		col = np.concatenate((col, range(size)))
150		data= np.concatenate((data,mult.flatten()))
151
152	nz = data!=0
153	return data[nz],(row[nz],col[nz])
154	
155	
156	
157	
158	
159