← Back to branch summary

~gabrib/cbc.solve/swingWithDolfin1.2

« back to all changes in this revision

Viewing changes to test/swing/fsinewton/unittests/test_jacobian.py

  • Committer: Gabriel Balaban
  • Date: 2013-08-05 02:30:48 UTC
  • Revision ID: gabrib@simula.no-20130805023048-5kbvvzay0q1dcpas
add

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/swing/fsinewton/unittests/test_jacobian.py
1
 
"""
2
 
A set of tests to insure that the manually computed jacobian is correct.
3
 
Testing is done against an automatically derived jacobian
4
 
using dolfin.derivative()
5
 
 
6
 
 
7
 
Fluid Block variables     U_F,P_F,L_U
8
 
Structure Block variables D_S,U_S
9
 
Mesh Block Variables      D_F,L_D
10
 
"""
11
 
 
12
 
__author__ = "Gabriel Balaban"
13
 
__copyright__ = "Copyright (C) 2012 Simula Research Laboratory and %s" % __author__
14
 
__license__  = "GNU GPL Version 3 or any later version"
15
 
 
16
 
import io
17
 
from dolfin import *
18
 
import numpy as np
19
 
import demo.swing.analytic.newtonanalytic as nana
20
 
from demo.swing.minimal.minimalproblem import FSIMini
21
 
import cbc.swing.fsinewton.solver.solver_fsinewton as sfsi
22
 
import cbc.swing.fsinewton.utils.misc_func as mf
23
 
import cbc.swing.fsinewton.solver.spaces as spaces
24
 
from cbc.swing.parameters import fsinewton_params
25
 
        
26
 
class TestJacobians(object):
27
 
    "Test to make sure that the manual jacobian matches the analytical one",
28
 
    def setup_class(self):
29
 
        self.problems = [nana.NewtonAnalytic()] #FSIMini()
30
 
 
31
 
        #Tolerance for diffreports
32
 
 
33
 
        #Fixme The tolerance should actually be 1.0e-14 but for some reason the the buffered jacobian
34
 
        #doesn't match at the higher tolerance level. From experience it seems to have good convergence
35
 
        #though.
36
 
        self.TOL = 1.0e-11
37
 
 
38
 
    def test_jacobians(self):
39
 
        "Test the manual FSI Jacobian against the analytic Jacobian"
40
 
        #Todo continue the test into a few Newton iterations        
41
 
        #Change the default solver params
42
 
        fsinewton_params["optimization"]["reuse_jacobian"] = False
43
 
        fsinewton_params["optimization"]["simplify_jacobian"] = False
44
 
        fsinewton_params["plot"] = False
45
 
        fsinewton_params["fluid_domain_time_discretization"] = "mid-point"
46
 
        
47
 
        for prob in self.problems:
48
 
            print "Testing Jacobian with problem", prob.__str__()
49
 
            fsisolvers = {}
50
 
            newtonsolvers = {}
51
 
            residuals = {}
52
 
            jacobians = {}
53
 
            blocks= {}
54
 
 
55
 
            #Check the jacobians and residuals pairwise
56
 
            pairs = [("auto","manual"),
57
 
                     ("auto","buff"),
58
 
                     ("manual","buff")]
59
 
 
60
 
            blocknames = ["J_FF","J_FS","J_FM","J_SF","J_SS","J_SM","J_MF","J_MS","J_MM"]
61
 
 
62
 
            t = "buff"
63
 
            fsinewton_params["solve"] = False
64
 
            fsisolvers[t] = sfsi.FSINewtonSolver(prob,fsinewton_params)
65
 
            fsisolvers[t].prepare_solve()
66
 
            newtonsolvers[t] = fsisolvers[t].newtonsolver
67
 
            newtonsolvers[t].build_residual()
68
 
            newtonsolvers[t].build_jacobian()
69
 
            
70
 
            for t1,t2 in pairs:
71
 
                info_blue("Testing Jacobians of %s %s \n \n"%(t1,t2))
72
 
            
73
 
                for t in t1,t2:
74
 
                    fsinewton_params["jacobian"] = t
75
 
                    #Create a fresh solver
76
 
                    fsisolvers[t] = sfsi.FSINewtonSolver(prob,fsinewton_params)
77
 
                    fsisolvers[t].prepare_solve()
78
 
                    newtonsolvers[t] = fsisolvers[t].newtonsolver
79
 
                    newtonsolvers[t].build_residual()
80
 
                    newtonsolvers[t].build_jacobian()
81
 
 
82
 
                #Do two newton iterations
83
 
                for i in range(2):
84
 
                    for t in t1,t2:
85
 
                        residuals[t] = newtonsolvers[t].F.array()
86
 
                        jacobians[t] = newtonsolvers[t].J.array()
87
 
                    
88
 
                        #Set up the Jacobian blocks, fsi space should always be the same
89
 
                        sublocator = spaces.FSISubSpaceLocator(fsisolvers[t1].spaces.fsispace)
90
 
 
91
 
                        #Get the blocks of the matrix
92
 
                        blocks[t] = self.fsiblocks(jacobians[t],sublocator)
93
 
                    
94
 
                    #Check that the residual is the same
95
 
                    print "Checking residuals"
96
 
                    diff = residuals[t1] - residuals[t2]
97
 
                    print np.all(diff < self.TOL)
98
 
                    assert np.all(diff < self.TOL),"Error in residuals %s %s"%(t1,t2)
99
 
                    
100
 
                    #Check each block
101
 
                    for i,blockname in enumerate(blocknames):
102
 
                        print "\nChecking block ",blockname
103
 
                        mf.diffmatrix_report(blocks[t1][i],blocks[t2][i],self.TOL)
104
 
                        diff = blocks[t1][i] - blocks[t2][i]
105
 
                        print  np.all(diff < self.TOL)
106
 
                        assert np.all(diff < self.TOL),\
107
 
                           "Error in jacobian '%s' '%s' comparison. Block %s doesn't match at TOL = %f"%(t1,t2,blockname,self.TOL)
108
 
                    newtonsolvers[t1].step(newtonsolvers[t1].tol)
109
 
                    newtonsolvers[t2].step(newtonsolvers[t2].tol)
110
 
                
111
 
    def fsiblocks(self,J,sl):
112
 
        """
113
 
        Divide the jacobian matrix into blocks according to the subspace locator sl
114
 
        Fluid Block variables     U_F,P_F,L_U
115
 
        Structure Block variables D_S,U_S
116
 
        Mesh Block Variables      D_F,L_D
117
 
        """
118
 
        sl.fluidend = sl.spaceends["L_U"]
119
 
        sl.strucend = sl.spaceends["U_S"]
120
 
        
121
 
        J_FF = J[:sl.fluidend,:sl.fluidend]
122
 
        J_FS = J[:sl.fluidend,sl.fluidend:sl.strucend]
123
 
        J_FM = J[:sl.fluidend,sl.strucend:]
124
 
        J_SF = J[sl.fluidend:sl.strucend,:sl.fluidend]
125
 
        J_SS = J[sl.fluidend:sl.strucend,sl.fluidend:sl.strucend]
126
 
        J_SM = J[sl.fluidend:sl.strucend,sl.strucend:]
127
 
        J_MF = J[sl.strucend:,:sl.fluidend]
128
 
        J_MS = J[sl.strucend:,sl.fluidend:sl.strucend]
129
 
        J_MM = J[sl.strucend:,sl.strucend:]
130
 
        return (J_FF,J_FS,J_FM,J_SF,J_SS,J_SM,J_MF,J_MS,J_MM)
131
 
 
132
 
if __name__ == "__main__":
133
 
    t = TestJacobians()
134
 
    t.setup_class()
135
 
    t.test_jacobians()