TL: first major rework of dedalus playground

Author: tlunet

pySDC/playgrounds/dedalus/__init__.py

@@ -0,0 +1 @@
+#

pySDC/playgrounds/dedalus/advection.py

@@ -0,0 +1,126 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Quick demo to solve the 1D advection-diffusion with Dedalus
+using the pySDC interface
+"""
+# Base user imports
+import numpy as np
+import matplotlib.pyplot as plt
+import dedalus.public as d3
+
+# pySDC imports
+from pySDC.playgrounds.dedalus.interface.problem import DedalusProblem
+from pySDC.playgrounds.dedalus.interface.sweeper import DedalusSweeperIMEX
+from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+
+# -----------------------------------------------------------------------------
+# User parameters
+# -----------------------------------------------------------------------------
+listK = [0, 1, 2]   # list of initial wavenumber in the solution (amplitude 1)
+nu = 1e-2           # viscosity (unitary velocity)
+
+# -- Space discretisation
+xEnd = 2*np.pi      # domain [0, xEnd]
+nX = 16             # number of points in x (periodic domain)
+
+# -- Time integration
+nSweeps = 4         # number of SDC iterations (sweeps)
+nNodes = 4          # number of SDC quadrature nodes
+tEnd = 2*np.pi      # final simulation time (starts at 0)
+nSteps = 50         # number of time-steps
+
+# -----------------------------------------------------------------------------
+# Solver setup
+# -----------------------------------------------------------------------------
+
+# -- Dedalus space grid
+coords = d3.CartesianCoordinates('x')
+dist = d3.Distributor(coords, dtype=np.float64)
+xbasis = d3.RealFourier(coords['x'], size=nX, bounds=(0, xEnd))
+u = dist.Field(name='u', bases=xbasis)
+
+# -- Initial solution
+x = xbasis.local_grid(dist, scale=1)
+u0 = np.sum([np.cos(k*x) for k in listK], axis=0)
+np.copyto(u['g'], u0)
+u0 = u.copy()   # store initial field for later
+
+# -- Problem
+dx = lambda f: d3.Differentiate(f, coords['x'])
+problem = d3.IVP([u], namespace=locals())
+problem.add_equation("dt(u) + dx(u) - nu*dx(dx(u)) = 0")
+
+nSweeps = 4
+nNodes = 4
+tEnd = 2*np.pi
+nSteps = 100
+dt = tEnd / nSteps
+
+# -- pySDC controller settings
+description = {
+    # Sweeper and its parameters
+    "sweeper_class": DedalusSweeperIMEX,
+    "sweeper_params": {
+        "quad_type": "RADAU-RIGHT",
+        "num_nodes": nNodes,
+        "node_type": "LEGENDRE",
+        "initial_guess": "copy",
+        "do_coll_update": False,
+        "QI": "IE",
+        "QE": "EE",
+        'skip_residual_computation':
+            ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE'),
+    },
+    # Step parameters
+    "step_params": {
+        "maxiter": 1,
+    },
+    # Level parameters
+    "level_params": {
+        "dt": dt,
+        "restol": -1,
+        "nsweeps": nSweeps,
+    },
+    "problem_class": DedalusProblem,
+    "problem_params": {
+        'problem': problem,
+        'nNodes': nNodes,
+    }
+}
+
+controller = controller_nonMPI(
+    num_procs=1, controller_params={'logger_level': 30},
+    description=description)
+
+prob = controller.MS[0].levels[0].prob
+uSol = prob.solver.state
+tVals = np.linspace(0, tEnd, nSteps + 1)
+
+for i in range(nSteps):
+    uSol, _ = controller.run(u0=uSol, t0=tVals[i], Tend=tVals[i + 1])
+
+# -- Plotting solution in real and Fourier space
+fig, (ax1, ax2) = plt.subplots(1, 2)
+fig.suptitle(rf"Advection-diffusion with $\nu={nu:1.1e}$")
+
+plt.sca(ax1)
+plt.title("Real space")
+plt.plot(u0['g'], label='$u_0$')
+plt.plot(u['g'], label='$u(T)$')
+plt.legend()
+plt.grid(True)
+plt.xlabel("$x$")
+plt.ylabel("$u$")
+
+plt.sca(ax2)
+plt.title("Coefficient space")
+plt.plot(u0['c'], 'o', mfc="none", label='$u_0$')
+plt.plot(u['c'], 's', mfc="none", ms=10, label='$u(t)$')
+plt.legend()
+plt.grid(True)
+plt.xlabel(r"$\kappa$")
+plt.ylabel("$u$")
+
+fig.set_size_inches(12, 5)
+plt.tight_layout()

pySDC/playgrounds/dedalus/demo.py

@@ -10,22 +10,24 @@
 import logging
 logger = logging.getLogger(__name__)
 
-from problem import DedalusProblem
-from sweeper import DedalusSweeperIMEX
+from pySDC.playgrounds.dedalus.interface.problem import DedalusProblem
+from pySDC.playgrounds.dedalus.interface.sweeper import DedalusSweeperIMEX
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 
-
-# Parameters
+# Space parameters
 Lx = 10
-Nx = 1024
+Nx = 512
 a = 1e-4
 b = 2e-4
 dealias = 3/2
 dtype = np.float64
 
+# Time-integration parameters
+nSweeps = 4
+nNodes = 4
 tEnd = 10
-nSteps = 10000
-
+nSteps = 5000
+timeStep = tEnd / nSteps
 
 # Bases
 xcoord = d3.Coordinate('x')
@@ -47,11 +49,6 @@
 n = 20
 u['g'] = np.log(1 + np.cosh(n)**2/np.cosh(n*(x-0.2*Lx))**2) / (2*n)
 
-# pySDC parameters
-dt = tEnd / nSteps
-nSweeps = 1
-nNodes = 4
-
 description = {
     # Sweeper and its parameters
     "sweeper_class": DedalusSweeperIMEX,
@@ -61,8 +58,8 @@
         "node_type": "LEGENDRE",
         "initial_guess": "copy",
         "do_coll_update": False,
-        "QI": "IE",
-        "QE": "EE",
+        "QI": "MIN-SR-S",
+        "QE": "PIC",
         'skip_residual_computation':
             ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE'),
     },
@@ -72,7 +69,7 @@
     },
     # Level parameters
     "level_params": {
-        "dt": dt,
+        "dt": timeStep,
         "restol": -1,
         "nsweeps": nSweeps,
     },
@@ -104,7 +101,6 @@
     if (i+1) % 100 == 0:
         print(f"step {i+1}/{nSteps}")
     if (i+1) % 25 == 0:
-
         u.change_scales(1)
         u_list.append(np.copy(u['g']))
         t_list.append(tVals[i])
@@ -121,4 +117,4 @@
 plt.ylabel('t')
 plt.title(f'KdV-Burgers, (a,b)=({a},{b})')
 plt.tight_layout()
-plt.savefig("KdV_Burgers_pySDC.pdf")
+plt.savefig("KdV_Burgers_interface.pdf")