вычисление скорости сходимости в фенике - PullRequest
Новая
       51

вычисление скорости сходимости в фенике

0 голосов
/ 13 ноября 2018

Я новичок в fenics, и мне действительно нужна помощь для отладки моего кода. По сути, я пытаюсь решить зависящее от времени уравнение в частных производных, и я хотел бы вычислить степень сходимости норм ошибок, которые у меня есть. однако я продолжаю получать эту ошибку, которую я не понимаю, что это значит или даже как ее преодолеть! 

rates[degree][error_type] = []

KeyError: 1

и это мой полный код 

T = 2.0       # final time
beta =  1.2   # parameter beta
tol = 3E-16
bound_a=0.0
bound_b=1.0
r=1.0
m=0.1


num_levels=3
max_degree=3
h = {}  # discretization parameter: h[degree][level]
E = {}  # error measure(s): E[degree][level][error_type]

# Iterate over degrees and mesh refinement levels
degrees = range(1, max_degree + 1)
for degree in degrees:
    h[degree] = []
    E[degree] = []
    for i in range(num_levels):
        num_steps = pow(4,i) # number of time steps
        dt = T / num_steps   # time step size
        # Create mesh and define function space
        mesh_division= int(sqrt(num_steps))
        mesh = IntervalMesh(mesh_division,bound_a,bound_b)
        h[degree].append(1.0 / mesh_division)
        V = FunctionSpace(mesh, 'P', 1)

        # Define boundary condition
        u_D = Expression('(b-a)/2 + beta*t +(x[0]-a)*(x[0]-b)', degree=1, beta=beta,a=bound_a,b=bound_b, t=0.0)

        def boundary(x, on_boundary):
            return on_boundary and (abs((x[0]-bound_a)*(x[0]-bound_b))<tol)


        bc = DirichletBC(V, u_D, boundary)

        # Define initial value
        u_n = interpolate(u_D, V)

        # Define variational problem
        u = TrialFunction(V)
        v = TestFunction(V)
        f = Constant(0.0)

        F = u*v*dx + dt*inner(grad(u), grad(v))*dx + dt*inner(grad(u_n), grad(u))*dx -dt*r*(m-u_n)*v*dx- (u_n + dt*f)*v*dx
        a, L = lhs(F), rhs(F)

        # Time-stepping
        u = Function(V)
        time = 0.0

        for n in range(num_steps):

            # Update current time
            time += dt
            u_D.t = time

            # Compute solution
            solve(a == L, u, bc)
            u_e= interpolate(u_D, V)

            def compute_errors(u_e, u):
                # Infinity norm based on nodal values
                E1 = abs(u_e.vector().get_local() - u.vector().get_local()).max()

                # L2 norm
                E2 = errornorm(u_e, u, norm_type='L2')

                # H1 seminorm
                E3 = errornorm(u_e, u, norm_type='H10')

                # Collect error measures in a dictionary with self-explanatory keys
                errors = {'infinity norm': E1,
                          'L2 norm': E2,
                          'H10 seminorm': E3}

                return errors

            errors = compute_errors(u_e, u)
            E[degree].append(errors)
            #Compute convergence rates
            etypes = list(E[1][0].keys())
            rates = {}
            for error_type in sorted(etypes):
                 rates[degree][error_type] = []
                 for i in range(1, num_levels):
                     Ei = E[degree][i][error_type]
                     r = ln(Ei / Eim1) / ln(h[degree][i] / h[degree][i - 1])
                     rates[degree][error_type].append(round(r, 2))

            # Update previous solution
            u_n.assign(u)

любая помощь будет оценена. Большое спасибо

...