Steady.py 3.19 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import numpy as np
import time
import scipy.sparse
import scipy.sparse.linalg

class Steady :
    def __init__(self, law, dof, jac, ilu=True) :
        self.law = law
        self.meshJacobian = jac
        self.dof = dof
        self.preconditioner = None
        self.tspsolve = 0
        self.tvt = 0
        self.tas = 0
        self.iterative = True
        self.ilu = ilu
        self.x = self.dof.new_vector()

    def solve(self, rtol = 1e-8, atol = 1e-10) :
        x = self.x
        self.dof.update_boundary_conditions(self.meshJacobian)
        x[self.dof._fixedidx] = self.dof._fixedval
        firstR = None
        newtonIter = 0
        ttot = time.clock()
        while True :
            tic = time.clock()
            res, jac = self.law.volumeTerm(self.meshJacobian, x)
            self.tvt += time.clock() - tic

            tic = time.clock()
            for row in self.dof._fixedidx :
                a = slice(jac.indptr[row], jac.indptr[row +1])
                jac.data[a] = np.where(jac.indices[a] == row, 1., 0.)
            jac.eliminate_zeros()
            self.tas += time.clock() - tic
            res[self.dof._fixedidx] = 0
            resn = np.linalg.norm(res)
            if firstR :
                print("newton (%i) % 6.2e %.1f" % (newtonIter, resn, -np.log10(resn/firstR)))
                if (resn/firstR < rtol or resn < atol) :
                    break
            else :
                firstR = resn
                print("newton (%i) % 6.2e" % (newtonIter, resn))
            newtonIter += 1
            if newtonIter > 10:
                break
            tic = time.clock()
            matrix = jac.tocsc()
            if not self.preconditioner :
                P = scipy.sparse.linalg.splu(matrix)
                self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
            if self.iterative :
                info = 1
                if self.preconditioner :
                    #I think there is a bug in scipy, without the callback argument, maxiter as no effect
                    dx, info = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=50, callback=lambda x:None)
                #if info != 0 and self.ilu:
                #    print("*** PRECOND ***")
                #    P = scipy.sparse.linalg.spilu(matrix)
                #    self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
                #    dx,info  = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=50, callback=lambda x:None)
                if info != 0 :
                    print("*** PRECOND FULL LU***")
                    P = scipy.sparse.linalg.splu(matrix)
                    self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
                    dx,info  = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=15, callback=lambda x:None)
            else :
                dx = scipy.sparse.linalg.spsolve(matrix, res)
            x -= dx
            self.tspsolve += time.clock() - tic

        print("spsolve: %.3g  term: %.3g  assemble: %.3g" % (self.tspsolve, self.tvt, self.tas))
        return x