Coverage for pySDC/tutorial/step_1/C_collocation_problem_setup.py: 100%

1import numpy as np 

2import scipy.sparse as sp 

3from pathlib import Path 

4 

5from pySDC.core.Collocation import CollBase 

6from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced 

7 

8 

9def main(): 

10 """ 

11 A simple test program to create and solve a collocation problem directly 

12 """ 

13 

14 # instantiate problem 

15 prob = heatNd_unforced( 

16 nvars=1023, # number of degrees of freedom 

17 nu=0.1, # diffusion coefficient 

18 freq=4, # frequency for the test value 

19 bc='dirichlet-zero', # boundary conditions 

20 ) 

21 

22 # instantiate collocation class, relative to the time interval [0,1] 

23 coll = CollBase(num_nodes=3, tleft=0, tright=1, node_type='LEGENDRE', quad_type='RADAU-RIGHT') 

24 

25 # set time-step size (warning: the collocation matrices are relative to [0,1], see above) 

26 dt = 0.1 

27 

28 # solve collocation problem 

29 err = solve_collocation_problem(prob=prob, coll=coll, dt=dt) 

30 

31 Path("data").mkdir(parents=True, exist_ok=True) 

32 f = open('data/step_1_C_out.txt', 'w') 

33 out = 'Error of the collocation problem: %8.6e' % err 

34 f.write(out + '\n') 

35 print(out) 

36 f.close() 

37 

38 assert err <= 4e-04, "ERROR: did not get collocation error as expected, got %s" % err 

39 

40 

41def solve_collocation_problem(prob, coll, dt): 

42 """ 

43 Routine to build and solve the linear collocation problem 

44 

45 Args: 

46 prob: a problem instance 

47 coll: a collocation instance 

48 dt: time-step size 

49 

50 Return: 

51 the analytic error of the solved collocation problem 

52 """ 

53 

54 # shrink collocation matrix: first line and column deals with initial value, not needed here 

55 Q = coll.Qmat[1:, 1:] 

56 

57 # build system matrix M of collocation problem 

58 M = sp.eye(prob.nvars[0] * coll.num_nodes) - dt * sp.kron(Q, prob.A) 

59 

60 # get initial value at t0 = 0 

61 u0 = prob.u_exact(t=0) 

62 # fill in u0-vector as right-hand side for the collocation problem 

63 u0_coll = np.kron(np.ones(coll.num_nodes), u0) 

64 # get exact solution at Tend = dt 

65 uend = prob.u_exact(t=dt) 

66 

67 # solve collocation problem directly 

68 u_coll = sp.linalg.spsolve(M, u0_coll) 

69 

70 # compute error 

71 err = np.linalg.norm(u_coll[-prob.nvars[0] :] - uend, np.inf) 

72 

73 return err 

74 

75 

76if __name__ == "__main__": 

77 main()