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

1from pathlib import Path 

2import numpy as np 

3import pickle 

4 

5from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

6from pySDC.projects.DAE.problems.simple_DAE import problematic_f 

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

8from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep 

9from pySDC.helpers.stats_helper import get_sorted 

10from pySDC.implementations.hooks.log_solution import LogSolution 

11 

12 

13def main(): 

14 """ 

15 A simple test program to see the fully implicit SDC solver in action 

16 """ 

17 # initialize level parameters 

18 level_params = dict() 

19 level_params['restol'] = 1e-6 

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

21 

22 # initialize sweeper parameters 

23 sweeper_params = dict() 

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

25 sweeper_params['num_nodes'] = 3 

26 

27 # initialize problem parameters 

28 problem_params = dict() 

29 problem_params['newton_tol'] = 1e-12 

30 

31 # initialize step parameters 

32 step_params = dict() 

33 step_params['maxiter'] = 40 

34 

35 # initialize controller parameters 

36 controller_params = dict() 

37 controller_params['logger_level'] = 30 

38 controller_params['hook_class'] = [LogSolution, LogGlobalErrorPostStep] 

39 

40 # Fill description dictionary for easy hierarchy creation 

41 description = dict() 

42 description['problem_class'] = problematic_f 

43 description['problem_params'] = problem_params 

44 description['sweeper_class'] = fully_implicit_DAE 

45 description['sweeper_params'] = sweeper_params 

46 description['level_params'] = level_params 

47 description['step_params'] = step_params 

48 

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

50 

51 # instantiate the controller 

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

53 

54 # set time parameters 

55 t0 = 0.0 

56 Tend = 1.0 

57 

58 # get initial values on finest level 

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

60 uinit = P.u_exact(t0) 

61 

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

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

64 

65 # check error 

66 err = get_sorted(stats, type='e_global_post_step', sortby='time') 

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

68 print(f"Error is {err}") 

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

70 

71 # store results 

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

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

74 sol_data = np.array([[sol[j][1][i] for j in range(len(sol))] for i in range(P.nvars)]) 

75 

76 data = dict() 

77 data['dt'] = sol_dt 

78 data['solution'] = sol_data 

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

80 

81 print("Done") 

82 

83 

84if __name__ == "__main__": 

85 main()