Coverage for pySDC / projects / Monodomain / problem_classes / TestODE.py: 91%

1import logging 

2import numpy as np 

3from pySDC.core.problem import Problem 

4from pySDC.core.common import RegisterParams 

5from pySDC.implementations.datatype_classes.mesh import mesh 

6from pySDC.projects.Monodomain.datatype_classes.my_mesh import imexexp_mesh 

7 

8""" 

9Here we define the problems classes for the multirate Dahlquist test equation y'=lambda_I*y + lambda_E*y + lambda_e*y 

10Things are done so that it is compatible witht the sweepers. 

11""" 

12 

13 

14class Parabolic(RegisterParams): 

15 def __init__(self, **problem_params): 

16 self._makeAttributeAndRegister(*problem_params.keys(), localVars=problem_params, readOnly=True) 

17 self.shape = (1,) 

18 self.init = ((1,), None, np.dtype("float64")) 

19 

20 

21class TestODE(Problem): 

22 def __init__(self, **problem_params): 

23 self.logger = logging.getLogger("step") 

24 

25 self.parabolic = Parabolic(**problem_params) 

26 self.size = 1 # one state variable 

27 self.init = ((self.size, *self.parabolic.init[0]), self.parabolic.init[1], self.parabolic.init[2]) 

28 

29 # invoke super init 

30 super(TestODE, self).__init__(self.init) 

31 # store all problem params dictionary values as attributes 

32 self._makeAttributeAndRegister(*problem_params.keys(), localVars=problem_params, readOnly=True) 

33 

34 # initial and end time 

35 self.t0 = 0.0 

36 self.Tend = 1.0 if self.end_time < 0.0 else self.end_time 

37 

38 # set lambdas, if not provided by user 

39 if not hasattr(self, 'lmbda_laplacian'): 

40 self.lmbda_laplacian = -5.0 

41 if not hasattr(self, 'lmbda_gating'): 

42 self.lmbda_gating = -10.0 

43 if not hasattr(self, 'lmbda_others'): 

44 self.lmbda_others = -1.0 

45 

46 self.dtype_u = mesh 

47 self.dtype_f = mesh 

48 

49 def init_exp_extruded(self, new_dim_shape): 

50 return ((*new_dim_shape, 1, self.init[0][1]), self.init[1], self.init[2]) 

51 

52 def initial_value(self): 

53 u0 = self.dtype_u(self.init, val=1.0) 

54 

55 return u0 

56 

57 def eval_f(self, u, t, fh=None): 

58 if fh is None: 

59 fh = self.dtype_f(init=self.init, val=0.0) 

60 

61 fh[0] = (self.lmbda_laplacian + self.lmbda_gating + self.lmbda_others) * u[0] 

62 

63 return fh 

64 

65 

66class MultiscaleTestODE(TestODE): 

67 def __init__(self, **problem_params): 

68 super(MultiscaleTestODE, self).__init__(**problem_params) 

69 

70 self.dtype_f = imexexp_mesh 

71 

72 self.rhs_stiff_indeces = [0] 

73 self.rhs_stiff_args = [0] 

74 self.rhs_nonstiff_indeces = [0] 

75 self.rhs_nonstiff_args = [0] 

76 self.rhs_exp_args = [0] 

77 self.rhs_exp_indeces = [0] 

78 self.rhs_non_exp_indeces = [] 

79 

80 self.constant_lambda_and_phi = True 

81 

82 self.one = self.dtype_u(init=self.init, val=1.0) 

83 

84 def solve_system(self, rhs, factor, u0, t, u_sol=None): 

85 if u_sol is None: 

86 u_sol = self.dtype_u(init=self.init, val=0.0) 

87 

88 u_sol[0] = rhs[0] / (1 - factor * self.lmbda_laplacian) 

89 

90 return u_sol 

91 

92 def eval_f(self, u, t, eval_impl=True, eval_expl=True, eval_exp=True, fh=None, zero_untouched_indeces=True): 

93 

94 if fh is None: 

95 fh = self.dtype_f(init=self.init, val=0.0) 

96 

97 if eval_expl: 

98 fh.expl[0] = self.lmbda_others * u[0] 

99 

100 if eval_impl: 

101 fh.impl[0] = self.lmbda_laplacian * u[0] 

102 

103 if eval_exp: 

104 fh.exp[0] = self.lmbda_gating * u[0] 

105 

106 return fh 

107 

108 def eval_lmbda_yinf_exp(self, u, lmbda, yinf): 

109 lmbda[0] = self.lmbda_gating 

110 yinf[0] = 0.0 

111 

112 def lmbda_eval(self, u, t, lmbda=None): 

113 if lmbda is None: 

114 lmbda = self.dtype_u(init=self.init, val=0.0) 

115 

116 lmbda[0] = self.lmbda_gating 

117 

118 return lmbda