Skip to content

Commit 801f1d2

Browse files
tlunettlunet
committed
TL: still refactoring
1 parent 7224da5 commit 801f1d2

File tree

6 files changed

+195
-82
lines changed

6 files changed

+195
-82
lines changed

‎pySDC/playgrounds/dedalus/demo_interface_burger.py‎

Lines changed: 11 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -3,24 +3,20 @@
33
"""
44
Demo script for the KdV-Burgers equation
55
"""
6-
76
import numpy as np
87
import matplotlib.pyplot as plt
9-
import dedalus.public as d3
108
import logging
119
logger = logging.getLogger(__name__)
1210

13-
from pySDC.playgrounds.dedalus.interface.problem import DedalusProblem
14-
from pySDC.playgrounds.dedalus.interface.sweeper import DedalusSweeperIMEX
11+
from pySDC.playgrounds.dedalus.problems import buildKdVBurgerProblem
12+
from pySDC.playgrounds.dedalus.interface import DedalusProblem, DedalusSweeperIMEX
1513
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
1614

1715
# Space parameters
18-
Lx = 10
19-
Nx = 512
20-
a = 1e-4
16+
xEnd = 10
17+
nX = 512
18+
nu = 1e-4
2119
b = 2e-4
22-
dealias = 3/2
23-
dtype = np.float64
2420

2521
# Time-integration parameters
2622
nSweeps = 4
@@ -29,25 +25,7 @@
2925
nSteps = 5000
3026
timeStep = tEnd / nSteps
3127

32-
# Bases
33-
xcoord = d3.Coordinate('x')
34-
dist = d3.Distributor(xcoord, dtype=dtype)
35-
xbasis = d3.RealFourier(xcoord, size=Nx, bounds=(0, Lx), dealias=dealias)
36-
37-
# Fields
38-
u = dist.Field(name='u', bases=xbasis)
39-
40-
# Substitutions
41-
dx = lambda A: d3.Differentiate(A, xcoord)
42-
43-
# Problem
44-
problem = d3.IVP([u], namespace=locals())
45-
problem.add_equation("dt(u) - a*dx(dx(u)) - b*dx(dx(dx(u))) = - u*dx(u)")
46-
47-
# Initial conditions
48-
x = dist.local_grid(xbasis)
49-
n = 20
50-
u['g'] = np.log(1 + np.cosh(n)**2/np.cosh(n*(x-0.2*Lx))**2) / (2*n)
28+
pData = buildKdVBurgerProblem(nX, xEnd, nu, b)
5129

5230
description = {
5331
# Sweeper and its parameters
@@ -75,12 +53,13 @@
7553
},
7654
"problem_class": DedalusProblem,
7755
"problem_params": {
78-
'problem': problem,
56+
'problem': pData["problem"],
7957
'nNodes': nNodes,
8058
}
8159
}
8260

8361
# Main loop
62+
u, x = [pData[key] for key in ["u", "x"]]
8463
u.change_scales(1)
8564
u_list = [np.copy(u['g'])]
8665
t_list = [0]
@@ -111,10 +90,10 @@
11190
plt.pcolormesh(
11291
x.ravel(), np.array(t_list), np.array(u_list), cmap='RdBu_r',
11392
shading='gouraud', rasterized=True, clim=(-0.8, 0.8))
114-
plt.xlim(0, Lx)
93+
plt.xlim(0, xEnd)
11594
plt.ylim(0, tEnd)
11695
plt.xlabel('x')
11796
plt.ylabel('t')
118-
plt.title(f'KdV-Burgers, (a,b)=({a},{b})')
97+
plt.title(r'KdV-Burgers, $(\nu,b)='f'({nu},{b})$')
11998
plt.tight_layout()
120-
plt.savefig("KdV_Burgers_interface.pdf")
99+
plt.savefig("demo_interface_burger.png")

‎pySDC/playgrounds/dedalus/demo_interface_advDiff.py‎ renamed to ‎pySDC/playgrounds/dedalus/demos/demo_interface_advDiff.py‎

Lines changed: 15 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -7,11 +7,10 @@
77
# Base user imports
88
import numpy as np
99
import matplotlib.pyplot as plt
10-
import dedalus.public as d3
1110

1211
# pySDC imports
13-
from pySDC.playgrounds.dedalus.interface.problem import DedalusProblem
14-
from pySDC.playgrounds.dedalus.interface.sweeper import DedalusSweeperIMEX
12+
from pySDC.playgrounds.dedalus.problems import buildAdvDiffProblem
13+
from pySDC.playgrounds.dedalus.interface import DedalusProblem, DedalusSweeperIMEX
1514
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
1615

1716
# -----------------------------------------------------------------------------
@@ -21,7 +20,6 @@
2120
nu = 1e-2 # viscosity (unitary velocity)
2221

2322
# -- Space discretisation
24-
xEnd = 2*np.pi # domain [0, xEnd]
2523
nX = 16 # number of points in x (periodic domain)
2624

2725
# -- Time integration
@@ -33,23 +31,7 @@
3331
# -----------------------------------------------------------------------------
3432
# Solver setup
3533
# -----------------------------------------------------------------------------
36-
37-
# -- Dedalus space grid
38-
coords = d3.CartesianCoordinates('x')
39-
dist = d3.Distributor(coords, dtype=np.float64)
40-
xbasis = d3.RealFourier(coords['x'], size=nX, bounds=(0, xEnd))
41-
u = dist.Field(name='u', bases=xbasis)
42-
43-
# -- Initial solution
44-
x = xbasis.local_grid(dist, scale=1)
45-
u0 = np.sum([np.cos(k*x) for k in listK], axis=0)
46-
np.copyto(u['g'], u0)
47-
u0 = u.copy() # store initial field for later
48-
49-
# -- Problem
50-
dx = lambda f: d3.Differentiate(f, coords['x'])
51-
problem = d3.IVP([u], namespace=locals())
52-
problem.add_equation("dt(u) + dx(u) - nu*dx(dx(u)) = 0")
34+
pData = buildAdvDiffProblem(nX, listK, nu)
5335

5436
nSweeps = 4
5537
nNodes = 4
@@ -84,7 +66,7 @@
8466
},
8567
"problem_class": DedalusProblem,
8668
"problem_params": {
87-
'problem': problem,
69+
'problem': pData["problem"],
8870
'nNodes': nNodes,
8971
}
9072
}
@@ -93,21 +75,28 @@
9375
num_procs=1, controller_params={'logger_level': 30},
9476
description=description)
9577

78+
# -----------------------------------------------------------------------------
79+
# Simulation run
80+
# -----------------------------------------------------------------------------
9681
prob = controller.MS[0].levels[0].prob
9782
uSol = prob.solver.state
9883
tVals = np.linspace(0, tEnd, nSteps + 1)
9984

10085
for i in range(nSteps):
10186
uSol, _ = controller.run(u0=uSol, t0=tVals[i], Tend=tVals[i + 1])
10287

103-
# -- Plotting solution in real and Fourier space
88+
# -----------------------------------------------------------------------------
89+
# Plotting solution in real and Fourier space
90+
# -----------------------------------------------------------------------------
91+
x, u, u0 = [pData[key] for key in ("x", "u", "u0")]
92+
10493
fig, (ax1, ax2) = plt.subplots(1, 2)
10594
fig.suptitle(rf"Advection-diffusion with $\nu={nu:1.1e}$")
10695

10796
plt.sca(ax1)
10897
plt.title("Real space")
109-
plt.plot(u0['g'], label='$u_0$')
110-
plt.plot(u['g'], label='$u(T)$')
98+
plt.plot(x, u0['g'], label='$u_0$')
99+
plt.plot(x, u['g'], label='$u(T)$')
111100
plt.legend()
112101
plt.grid(True)
113102
plt.xlabel("$x$")
@@ -124,3 +113,4 @@
124113

125114
fig.set_size_inches(12, 5)
126115
plt.tight_layout()
116+
plt.savefig("demo_interface_advDiff.png")

‎pySDC/playgrounds/dedalus/demo_timestepper_advection.py‎ renamed to ‎pySDC/playgrounds/dedalus/demos/demo_timestepper_advection.py‎

Lines changed: 14 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -8,51 +8,38 @@
88
import matplotlib.pyplot as plt
99
import dedalus.public as d3
1010

11+
from pySDC.playgrounds.dedalus.problems import buildAdvDiffProblem
1112
from pySDC.playgrounds.dedalus.timestepper import SpectralDeferredCorrectionIMEX
1213

1314

1415
# -----------------------------------------------------------------------------
1516
# User parameters
1617
# -----------------------------------------------------------------------------
1718
listK = [0, 1, 2] # list of initial wavenumber in the solution (amplitude 1)
18-
19-
# -- Space discretisation
20-
xEnd = 2*np.pi # domain [0, xEnd]
21-
nX = 64 # number of points in x (periodic domain)
22-
19+
nX = 128 # number of points in x (periodic domain)
2320

2421
# -----------------------------------------------------------------------------
2522
# Solver setup
2623
# -----------------------------------------------------------------------------
27-
28-
# -- Bases and field
29-
coords = d3.CartesianCoordinates('x')
30-
dist = d3.Distributor(coords, dtype=np.float64)
31-
xbasis = d3.RealFourier(coords['x'], size=nX, bounds=(0, xEnd))
32-
u = dist.Field(name='u', bases=xbasis)
33-
34-
# -- Initial solution
35-
x = xbasis.local_grid(dist, scale=1)
36-
u0 = np.sum([np.cos(k*x) for k in listK], axis=0)
37-
np.copyto(u['g'], u0)
38-
uInit = u.copy() # copy initial field for latter
39-
40-
# -- Problem
41-
dx = lambda f: d3.Differentiate(f, coords['x'])
42-
problem = d3.IVP([u], namespace=locals())
43-
problem.add_equation("dt(u) + dx(u) = 0")
24+
pData = buildAdvDiffProblem(nX, listK)
25+
u = pData["u"]
26+
uInit = pData["u0"]
27+
u0 = uInit['g'].data
4428

4529
# -- Prepare plots
4630
orderPlot = {'RK111': 1,
4731
'RK222': 2,
4832
'RK443': 3,
4933
'ERK4': 4}
5034

51-
plt.figure('Error')
5235
SpectralDeferredCorrectionIMEX.setParameters(
53-
nNodes=3, nodeType='LEGENDRE', quadType='RADAU-RIGHT',
36+
nNodes=4, nodeType='LEGENDRE', quadType='RADAU-RIGHT',
5437
implSweep='MIN-SR-FLEX', explSweep='PIC', initSweep='COPY')
5538

39+
# -----------------------------------------------------------------------------
40+
# Simulations
41+
# -----------------------------------------------------------------------------
42+
plt.figure('Error')
5643
for timeStepper in [d3.RK111, d3.RK222, d3.RK443, 1, 2, 3]:
5744

5845
useSDC = False
@@ -71,6 +58,7 @@
7158
print(f" -- solving using {scheme}")
7259

7360
# -- Build solver
61+
problem = pData["problem"]
7462
solver = problem.build_solver(timeStepper)
7563
solver.stop_sim_time = 2*np.pi
7664
name = timeStepper.__name__
@@ -106,6 +94,7 @@ def getErr(nStep):
10694
plt.legend()
10795
plt.grid(True)
10896
plt.tight_layout()
97+
plt.savefig("demo_timestepper_advDiff_convergence.png")
10998

11099
# Plotting solution in real and Fourier space
111100
fig, (ax1, ax2) = plt.subplots(1, 2)
@@ -131,3 +120,4 @@ def getErr(nStep):
131120

132121
fig.set_size_inches(12, 5)
133122
plt.tight_layout()
123+
plt.savefig("demo_timestepper_advDiff.png")

‎pySDC/playgrounds/dedalus/interface/__init__.py‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,7 @@
1-
#
1+
from pySDC.playgrounds.dedalus.interface.problem import DedalusProblem
2+
from pySDC.playgrounds.dedalus.interface.sweeper import DedalusSweeperIMEX
3+
4+
__all__ = [
5+
"DedalusProblem",
6+
"DedalusSweeperIMEX"
7+
]

‎pySDC/playgrounds/dedalus/problems/__init__.py‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
from pySDC.playgrounds.dedalus.problems.periodic1D import buildAdvDiffProblem, buildKdVBurgerProblem
2+
3+
__all__ = [
4+
"buildAdvDiffProblem",
5+
"buildKdVBurgerProblem"
6+
]

0 commit comments

Comments
 (0)