Coverage for pySDC/projects/DAE/run/synchronous_machine_playground.py: 100%

1from pathlib import Path 

2import numpy as np 

3import pickle 

4import statistics 

5 

6from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

7from pySDC.projects.DAE.problems.synchronous_machine import synchronous_machine_infinite_bus 

8from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE 

9from pySDC.projects.DAE.misc.HookClass_DAE import LogGlobalErrorPostStepDifferentialVariable 

10from pySDC.helpers.stats_helper import get_sorted 

11from pySDC.helpers.stats_helper import filter_stats 

12from pySDC.implementations.hooks.log_solution import LogSolution 

13 

14 

15def main(): 

16 """ 

17 A testing ground for the synchronous machine model 

18 """ 

19 # initialize level parameters 

20 level_params = dict() 

21 level_params['restol'] = 1e-7 

22 level_params['dt'] = 1e-1 

23 

24 # initialize sweeper parameters 

25 sweeper_params = dict() 

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

27 sweeper_params['num_nodes'] = 3 

28 sweeper_params['QI'] = 'LU' 

29 

30 # initialize problem parameters 

31 problem_params = dict() 

32 problem_params['newton_tol'] = 1e-3 

33 

34 # initialize step parameters 

35 step_params = dict() 

36 step_params['maxiter'] = 100 

37 

38 # initialize controller parameters 

39 controller_params = dict() 

40 controller_params['logger_level'] = 30 

41 controller_params['hook_class'] = [LogGlobalErrorPostStepDifferentialVariable, LogSolution] 

42 

43 # Fill description dictionary for easy hierarchy creation 

44 description = dict() 

45 description['problem_class'] = synchronous_machine_infinite_bus 

46 description['problem_params'] = problem_params 

47 description['sweeper_class'] = fully_implicit_DAE 

48 description['sweeper_params'] = sweeper_params 

49 description['level_params'] = level_params 

50 description['step_params'] = step_params 

51 

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

53 

54 # instantiate the controller 

55 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) 

56 

57 # set time parameters 

58 t0 = 0.0 

59 Tend = 1.0 

60 

61 # get initial values on finest level 

62 P = controller.MS[0].levels[0].prob 

63 uinit = P.u_exact(t0) 

64 

65 # call main function to get things done... 

66 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) 

67 

68 # check error (only available if reference solution was provided) 

69 # err = get_sorted(stats, type='e_global_differential_post_step', sortby='time') 

70 # err = np.linalg.norm([err[i][1] for i in range(len(err))], np.inf) 

71 # print(f"Error is {err}") 

72 

73 uend_ref = P.dtype_u(P.init) 

74 uend_ref.diff[:8] = ( 

75 8.30823565e-01, 

76 -4.02584174e-01, 

77 1.16966755e00, 

78 9.47592808e-01, 

79 -3.68076863e-01, 

80 -3.87492326e-01, 

81 3.10281509e-01, 

82 9.94039645e-01, 

83 ) 

84 uend_ref.alg[:6] = ( 

85 -7.77837831e-01, 

86 -1.67347611e-01, 

87 1.34810867e00, 

88 5.46223705e-04, 

89 1.29690691e-02, 

90 -8.00823474e-02, 

91 ) 

92 err = abs(uend.diff - uend_ref.diff) 

93 assert np.isclose(err, 0, atol=1e-4), "Error too large." 

94 

95 # store results 

96 sol = get_sorted(stats, type='u', sortby='time') 

97 sol_dt = np.array([sol[i][0] for i in range(len(sol))]) 

98 sol_data = np.array( 

99 [ 

100 [(sol[j][1].diff[id], sol[j][1].alg[ia]) for j in range(len(sol))] 

101 for id, ia in zip(range(len(uend.diff)), range(len(uend.alg))) 

102 ] 

103 ) 

104 niter = filter_stats(stats, type='niter') 

105 niter = np.fromiter(niter.values(), int) 

106 

107 data = dict() 

108 data['dt'] = sol_dt 

109 data['solution'] = sol_data 

110 data['niter'] = round(statistics.mean(niter)) 

111 pickle.dump(data, open("data/dae_conv_data.p", 'wb')) 

112 

113 print("Done") 

114 

115 

116if __name__ == "__main__": 

117 main()