scipylsys.py 3.52 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
4
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
6
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
7
8
9
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
12
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
13
14
15
16
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
17
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
18
# You should have received a copy of the GNU Lesser General Public License
19
# along with this program (see COPYING and COPYING.LESSER files).  If not,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
# see <http://www.gnu.org/licenses/>.
21

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
22
23
24
import numpy as np
import scipy.sparse
import scipy.sparse.linalg
25
from ._tools import timeit
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
26
27

class LinearSystem :
28
29

    def __init__(self, elements, n_fields,options, constraints = []) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
30
31
32
33
34
35
36
37
38
        localsize = n_fields*elements.shape[1]
        self.el = elements
        idx = (elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1,localsize])
        self.rows = np.repeat(idx[:,:,None],localsize,axis=2)
        self.rows = self.rows.reshape([-1])
        self.cols = np.repeat(idx[:,None,:],localsize,axis=1)
        self.cols = self.cols.reshape([-1])
        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
        self.localsize = localsize
39
40
41
42
43
44
45
46
47
48
        self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
        self.ndof = (np.max(elements) + 1)*n_fields
        self.globalsize = self.ndof
        for i, constraint in enumerate(self.constraints):
            self.rows = np.hstack((self.rows,
                                      constraint,
                                      np.full((constraint.shape[0],),i+self.ndof)))
            self.cols = np.hstack((self.cols,
                                      np.full((constraint.shape[0],),i+self.ndof),
                                      constraint))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
49

50
51
52
53
54
55
56
57
58
    @timeit
    def local_to_global(self,localv,localm, u, constraints_value = []):
        rhs = np.zeros((self.ndof + len(constraints_value,)))
        np.add.at(rhs,self.idx,localv)
        localm = localm.reshape([-1])
        for i, (c, cv) in enumerate(zip(self.constraints, constraints_value)):
            localm = np.hstack((localm, cv[0], cv[0]))
            rhs[i+self.ndof] = cv[1] + np.sum(u[c]*cv[0])
        coo = scipy.sparse.coo_matrix((localm, (self.rows,self.cols)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
59
        self.matrix = coo.tocsc()
60
61
62
63
64
65
        return rhs

    def solvelu(self,rhs) :
        scipy.sparse.linalg.use_solver(useUmfpack=False,assumeSortedIndices=True)
        r =  scipy.sparse.linalg.spsolve(self.matrix,rhs,permc_spec="MMD_ATA")[:self.ndof]
        return r
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
66

67
68
69
70
71
72
73
74
75
76
    def solvegmres(self,rhs) :
        print("ilu")
        #ilu = scipy.sparse.linalg.spilu(self.matrix,fill_factor=0)
        ilu = scipy.sparse.linalg.spilu(self.matrix,drop_tol=1e-5)
        print("ilu solve pc")
        pc = scipy.sparse.linalg.LinearOperator(self.matrix.shape,ilu.solve)
        print("gmres")
        return scipy.sparse.linalg.gmres(self.matrix,rhs,M=pc,restart=30,maxiter=1000)[0][:self.ndof]
        #return scipy.sparse.linalg.splu(self.matrix).solve(rhs)[:self.ndof]
        #return r
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
77

78
79
80
    @timeit
    def solve(self,rhs):
        return self.solvelu(rhs)