diffusion_project/part_d.py at master · gba/diffusion_project · GitHub

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
"""
The non-linear diffusion equation

rho*du/dt = \nabla \dot (alpha(u)\nabla u) + f(x, t),

with x \in Omega, t \in (0, T], is solved with only
one Picard iteration at every time step.

The PDE is discretized in time using the Backward
Euler method.

In this file we test the convergence rate of a
simple linear problem.

"""
import scitools.std as sci
import nose.tools as nt
from dolfin import *
import numpy, sys

Ct = 1.0
Cx = 0.01

h_values = [0.05, 0.01, 0.005, 0.001]
errors = []

file_1 = open('x_u_d3.dat','w')
file_e = open('x_ue_d3.dat','w')

for h in h_values:

    #     degree:        Degree of finite element
    #                    approximation
    #     n_elements:    List of number of elements
    #                    for each dimension
    #     dim:           Physical space dimension
    #     ref_domain:    Type of reference domain
    #     V:             Test function space
    degree = 1
    Nx = int(round(1.0/numpy.sqrt(Cx*h)))
    n_elements = [Nx, Nx]
    dim = len(n_elements)
    ref_domain = [UnitInterval, UnitSquare, UnitCube]
    mesh = ref_domain[dim-1](*n_elements)
    V = FunctionSpace(mesh, 'Lagrange', degree)

    #     The Neumann boundary conditions are
    #     implicitly included in the variational
    #     formulation.
    #
    #     Initial condition
    u0 = Expression('exp(-pi*pi*t)*cos(pi*x[0])', t=0.0)
    u_1 = interpolate(u0, V)

    #     dt:       Time step length
    #     T:        Stopping time for simulation
    T = 0.05
    Nt = int(round(T/(Ct*h)))
    dt = T/Nt
    print 'Nx = ', Nx, ', dt = ', dt

    #     Diffusion coefficient
    def alpha(u):
        return 1.0

    #     Source function
    f = Constant(0.0)
    #     Constant rho
    rho = 1.0

    #     Set up the variational formulation
    u = TrialFunction(V)
    v = TestFunction(V)
    a = (rho*inner(u, v) \
             + dt*inner(alpha(u_1)*nabla_grad(u), \
                            nabla_grad(v)))*dx
    L = (dt*inner(f, v) + rho*inner(u_1, v))*dx

    u = Function(V)
    t = dt
    #     Loop over time
    while t <= T:
        f.t = t
        #     Assemble the matrix and vector
        A = assemble(a)
        b = assemble(L)
        #     Solve the linear system of equations
        solve(A, u.vector(), b)

        #     Test the result against the exact solution
        u0.t = t
        u_e = interpolate(u0, V)
        e = u_e.vector().array() - u.vector().array()
        E = numpy.sqrt(numpy.sum(e**2) \
                           /u.vector().array().size)

        if t+dt > T:
            errors.append(E)
            #plot(u)
            #interactive()

            #     Write numerical and exact solutions
            #     to file
            # uv = u.vector().array()
            # x = numpy.zeros(len(uv))
            # for i in range(len(uv)):
            #     file_1.write(str(mesh.coordinates()[i][0])
            #                  +' '+str(u.vector().array()[i])+'\n')
            #     file_e.write(str(mesh.coordinates()[i][0])
            #                  +' '+str(u_e.vector().array()[i])+'\n')

        t += dt
        u_1.assign(u)

print 'errors = ', errors
nh = len(h_values)

#     Calculate convergence rates
r = [numpy.log(errors[i-1]/errors[i])/ \
         numpy.log(h_values[i-1]/h_values[i]) \
         for i in range(1, nh, 1)]
#     Print convergence parameter h and convergence
#     rate r
for i in range(1, nh):
    print h_values[i-1], r[i-1]

#     Check that convergence rate is 1 with a nose
#     test
diff = abs(r[nh-2]-1.0)
print 'diff = ', diff
nt.assert_almost_equal(diff, 0, delta=1E-2)

file_1.close()
file_e.close()