|
| 1 | +import numpy as np |
| 2 | +from pathlib import Path |
| 3 | +import matplotlib.pyplot as plt |
| 4 | + |
| 5 | +from pySDC.helpers.stats_helper import get_sorted |
| 6 | +from pySDC.projects.DAE.sweepers.SemiImplicitDAE import SemiImplicitDAE |
| 7 | +from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE |
| 8 | +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI |
| 9 | + |
| 10 | +from pySDC.implementations.hooks.log_solution import LogSolution |
| 11 | + |
| 12 | + |
| 13 | +def simulateDAE(): |
| 14 | + r""" |
| 15 | + Main function where things will be done. Here, the problem class ``DiscontinuousTestDAE`` is simulated using |
| 16 | + the ``SemiImplicitDAE`` sweeper, where only the differential variable is integrated using spectral quadrature. |
| 17 | + The problem usually contains a discrete event, but the simulation interval is chosen so that the state function |
| 18 | + does not have a sign change yet (the discrete event does not yet occur). Thus, the solution is smooth. |
| 19 | + """ |
| 20 | + # initialize level parameters |
| 21 | + level_params = { |
| 22 | + 'restol': 1e-12, |
| 23 | + 'dt': 0.1, |
| 24 | + } |
| 25 | + |
| 26 | + # initialize problem parameters |
| 27 | + problem_params = { |
| 28 | + 'newton_tol': 1e-6, |
| 29 | + } |
| 30 | + |
| 31 | + # initialize sweeper parameters |
| 32 | + sweeper_params = { |
| 33 | + 'quad_type': 'RADAU-RIGHT', |
| 34 | + 'num_nodes': 3, |
| 35 | + 'QI': 'IE', |
| 36 | + 'initial_guess': 'spread', |
| 37 | + } |
| 38 | + |
| 39 | + # initialize step parameters |
| 40 | + step_params = { |
| 41 | + 'maxiter': 60, |
| 42 | + } |
| 43 | + |
| 44 | + # initialize controller parameters |
| 45 | + controller_params = { |
| 46 | + 'logger_level': 30, |
| 47 | + 'hook_class': [LogSolution], |
| 48 | + } |
| 49 | + |
| 50 | + # fill description dictionary for easy step instantiation |
| 51 | + description = { |
| 52 | + 'problem_class': DiscontinuousTestDAE, |
| 53 | + 'problem_params': problem_params, |
| 54 | + 'sweeper_class': SemiImplicitDAE, |
| 55 | + 'sweeper_params': sweeper_params, |
| 56 | + 'level_params': level_params, |
| 57 | + 'step_params': step_params, |
| 58 | + } |
| 59 | + |
| 60 | + # instantiate controller |
| 61 | + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) |
| 62 | + |
| 63 | + t0 = 1.0 |
| 64 | + Tend = 3.0 |
| 65 | + |
| 66 | + P = controller.MS[0].levels[0].prob |
| 67 | + uinit = P.u_exact(t0) |
| 68 | + |
| 69 | + Path("data").mkdir(parents=True, exist_ok=True) |
| 70 | + |
| 71 | + # call main function to get things done... |
| 72 | + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) |
| 73 | + |
| 74 | + plotSolution(stats) |
| 75 | + |
| 76 | + |
| 77 | +def plotSolution(stats): |
| 78 | + r""" |
| 79 | + Here, the solution of the DAE is plotted. |
| 80 | +
|
| 81 | + Parameters |
| 82 | + ---------- |
| 83 | + stats : dict |
| 84 | + Contains all the statistics of the simulation. |
| 85 | + """ |
| 86 | + |
| 87 | + u_val = get_sorted(stats, type='u', sortby='time') |
| 88 | + t = np.array([me[0] for me in u_val]) |
| 89 | + y = np.array([me[1].diff[0] for me in u_val]) |
| 90 | + z = np.array([me[1].alg[0] for me in u_val]) |
| 91 | + |
| 92 | + plt.figure(figsize=(8.5, 8.5)) |
| 93 | + plt.plot(t, y, label='Differential variable y') |
| 94 | + plt.plot(t, z, label='Algebraic variable z') |
| 95 | + plt.legend(frameon=False, fontsize=12, loc='upper left') |
| 96 | + plt.xlabel(r"Time $t$") |
| 97 | + plt.ylabel(r"Solution $y(t)$, $z(t)$") |
| 98 | + |
| 99 | + plt.savefig('data/solution.png', dpi=300, bbox_inches='tight') |
| 100 | + plt.close() |
| 101 | + |
| 102 | + |
| 103 | +if __name__ == "__main__": |
| 104 | + simulateDAE() |
0 commit comments