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

21 statements  

« prev     ^ index     » next       coverage.py v7.6.9, created at 2024-12-20 14:51 +0000

1import numpy as np 

2from pathlib import Path 

3 

4from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced 

5 

6 

7def main(): 

8 """ 

9 A simple test program to set up a spatial problem and play with it 

10 """ 

11 # instantiate problem 

12 prob = heatNd_unforced( 

13 nvars=1023, # number of degrees of freedom 

14 nu=0.1, # diffusion coefficient 

15 freq=4, # frequency for the test value 

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

17 ) 

18 

19 # run accuracy test, get error back 

20 err = run_accuracy_check(prob) 

21 

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

23 f = open('data/step_1_A_out.txt', 'w') 

24 out = 'Error of the spatial accuracy test: %8.6e' % err 

25 f.write(out) 

26 print(out) 

27 f.close() 

28 

29 assert err <= 2e-04, "ERROR: the spatial accuracy is higher than expected, got %s" % err 

30 

31 

32def run_accuracy_check(prob): 

33 """ 

34 Routine to check the error of the Laplacian vs. its FD discretization 

35 

36 Args: 

37 prob: a problem instance 

38 

39 Returns: 

40 the error between the analytic Laplacian and the computed one of a given function 

41 """ 

42 

43 # create x values, use only inner points 

44 xvalues = np.array([(i + 1) * prob.dx for i in range(prob.nvars[0])]) 

45 

46 # create a mesh instance and fill it with a sine wave 

47 u = prob.dtype_u(init=prob.init) 

48 u[:] = np.sin(np.pi * prob.freq[0] * xvalues) 

49 

50 # create a mesh instance and fill it with the Laplacian of the sine wave 

51 u_lap = prob.dtype_u(init=prob.init) 

52 u_lap[:] = -((np.pi * prob.freq[0]) ** 2) * prob.nu * np.sin(np.pi * prob.freq[0] * xvalues) 

53 

54 # compare analytic and computed solution using the eval_f routine of the problem class 

55 err = abs(prob.eval_f(u, 0) - u_lap) 

56 

57 return err 

58 

59 

60if __name__ == "__main__": 

61 main()