Coverage for pySDC / projects / StroemungsRaum / sweepers / imex_1st_order_mass_NSE.py: 92%

1from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass 

2 

3 

4class imex_1st_order_mass_NSE(imex_1st_order_mass): 

5 """ 

6 This sweeper extends imex_1st_order_mass for a projection-based discretization of the incompressible 

7 Navier–Stokes equations. This requires passing the zero-to-node step size at each collocation 

8 node and handling the update of the velocity and the pressure within `solve_system` 

9 

10 """ 

11 

12 def update_nodes(self): 

13 """ 

14 Update the u- and f-values at the collocation nodes -> corresponds to a single sweep over all nodes 

15 

16 Returns: 

17 None 

18 """ 

19 # get current level and problem description 

20 L = self.level 

21 P = L.prob 

22 

23 # only if the level has been touched before 

24 assert L.status.unlocked 

25 

26 # get number of collocation nodes for easier access 

27 M = self.coll.num_nodes 

28 

29 # gather all terms which are known already (e.g. from the previous iteration) 

30 # this corresponds to u0 + QF(u^k) - QIFI(u^k) - QEFE(u^k) + tau 

31 

32 # get QF(u^k) 

33 integral = self.integrate() 

34 

35 # This is somewhat ugly, but we have to apply the mass matrix on u0 only on the finest level 

36 if L.level_index == 0: 

37 u0 = P.apply_mass_matrix(L.u[0]) 

38 else: 

39 u0 = L.u[0] 

40 

41 for m in range(M): 

42 # subtract QIFI(u^k)_m + QEFE(u^k)_m 

43 for j in range(M + 1): 

44 integral[m] -= L.dt * (self.QI[m + 1, j] * L.f[j].impl + self.QE[m + 1, j] * L.f[j].expl) 

45 # add initial value 

46 integral[m] += u0 

47 # add tau if associated 

48 if L.tau[m] is not None: 

49 integral[m] += L.tau[m] 

50 

51 # do the sweep 

52 for m in range(0, M): 

53 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

54 rhs = P.dtype_u(integral[m]) 

55 for j in range(m + 1): 

56 rhs += L.dt * (self.QI[m + 1, j] * L.f[j].impl + self.QE[m + 1, j] * L.f[j].expl) 

57 

58 L.u[m + 1], P.pn = P.solve_system( 

59 rhs, 

60 L.dt * self.QI[m + 1, m + 1], 

61 L.u[m + 1], 

62 L.time + L.dt * self.coll.nodes[m], 

63 L.dt * self.coll.nodes[m], 

64 ) 

65 

66 # update function values 

67 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) 

68 

69 # indicate presence of new values at this level 

70 L.status.updated = True 

71 return None