Coverage for pySDC/tutorial/step_2/B_my_first_sweeper.py: 100%

1from pathlib import Path 

2 

3from pySDC.core.step import Step 

4 

5from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced 

6from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

7 

8 

9def main(): 

10 """ 

11 A simple test program to run IMEX SDC for a single time step 

12 """ 

13 # initialize level parameters 

14 level_params = dict() 

15 level_params['restol'] = 1e-10 

16 level_params['dt'] = 0.1 

17 

18 # initialize sweeper parameters 

19 sweeper_params = dict() 

20 sweeper_params['quad_type'] = 'RADAU-RIGHT' 

21 sweeper_params['num_nodes'] = 3 

22 

23 # initialize problem parameters 

24 problem_params = dict() 

25 problem_params['nu'] = 0.1 # diffusion coefficient 

26 problem_params['freq'] = 4 # frequency for the test value 

27 problem_params['nvars'] = 1023 # number of degrees of freedom 

28 problem_params['bc'] = 'dirichlet-zero' # boundary conditions 

29 

30 # initialize step parameters 

31 step_params = dict() 

32 step_params['maxiter'] = 20 

33 

34 # Fill description dictionary for easy hierarchy creation 

35 description = dict() 

36 description['problem_class'] = heatNd_forced 

37 description['problem_params'] = problem_params 

38 description['sweeper_class'] = imex_1st_order 

39 description['sweeper_params'] = sweeper_params 

40 description['level_params'] = level_params 

41 description['step_params'] = step_params 

42 

43 # instantiate the step we are going to work on 

44 S = Step(description=description) 

45 

46 # run IMEX SDC test and check error, residual and number of iterations 

47 err, res, niter = run_imex_sdc(S) 

48 print('Error and residual: %12.8e -- %12.8e' % (err, res)) 

49 

50 assert err <= 1e-5, "ERROR: IMEX SDC iteration did not reduce the error enough, got %s" % err 

51 assert res <= level_params['restol'], "ERROR: IMEX SDC iteration did not reduce the residual enough, got %s" % res 

52 assert niter <= 12, "ERROR: IMEX SDC took too many iterations, got %s" % niter 

53 

54 

55def run_imex_sdc(S): 

56 """ 

57 Routine to run IMEX SDC on a single time step 

58 

59 Args: 

60 S: an instance of a step representing the time step 

61 

62 Returns: 

63 the error of SDC vs. exact solution 

64 the residual after the SDC sweeps 

65 the number of iterations 

66 """ 

67 # make shortcuts for the level and the problem 

68 L = S.levels[0] 

69 P = L.prob 

70 

71 # set initial time in the status of the level 

72 L.status.time = 0.1 

73 # compute initial value (using the exact function here) 

74 L.u[0] = P.u_exact(L.time) 

75 

76 # access the sweeper's predict routine to get things started 

77 # if we don't do this, the values at the nodes are not initialized 

78 L.sweep.predict() 

79 # compute the residual (we may be done already!) 

80 L.sweep.compute_residual() 

81 

82 # reset iteration counter 

83 S.status.iter = 0 

84 # run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough 

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

86 f = open('data/step_2_B_out.txt', 'w') 

87 while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol: 

88 # this is where the nodes are actually updated according to the SDC formulas 

89 L.sweep.update_nodes() 

90 # compute/update the residual 

91 L.sweep.compute_residual() 

92 # increment the iteration counter 

93 S.status.iter += 1 

94 out = 'Time %4.2f of %s -- Iteration: %2i -- Residual: %12.8e' % ( 

95 L.time, 

96 L.level_index, 

97 S.status.iter, 

98 L.status.residual, 

99 ) 

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

101 print(out) 

102 f.close() 

103 

104 # compute the interval's endpoint: this (and only this) will set uend, depending on the collocation nodes 

105 L.sweep.compute_end_point() 

106 # update the simulation time 

107 L.status.time += L.dt 

108 

109 # compute exact solution and compare 

110 uex = P.u_exact(L.status.time) 

111 err = abs(uex - L.uend) 

112 

113 return err, L.status.residual, S.status.iter 

114 

115 

116if __name__ == "__main__": 

117 main()