Skip to content

Commit

Permalink
Implemented changes requested by @tlunet
Browse files Browse the repository at this point in the history
  • Loading branch information
@brownbaerchen
brownbaerchen committed Sep 15, 2024
1 parent f33a7e9 commit 0e8e997
Show file tree
Hide file tree
Showing 9 changed files with 610 additions and 161 deletions.
435 changes: 384 additions & 51 deletions pySDC/helpers/spectral_helper.py
View file

Large diffs are not rendered by default.

34 changes: 28 additions & 6 deletions pySDC/implementations/problem_classes/Burgers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,18 @@ class Burgers1D(GenericSpectralLinear):
dtype_f = imex_mesh

def __init__(self, N=64, epsilon=0.1, BCl=1, BCr=-1, f=0, mode='T2U', **kwargs):
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
N (int): Spatial resolution
epsilon (float): viscosity
BCl (float): Value at left boundary
BCr (float): Value at right boundary
f (int): Frequency of the initial conditions
mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
"""
self._makeAttributeAndRegister('N', 'epsilon', 'BCl', 'BCr', 'f', 'mode', localVars=locals(), readOnly=True)

bases = [{'base': 'cheby', 'N': N}]
components = ['u', 'ux']
Expand Down Expand Up @@ -96,7 +107,7 @@ def eval_f(self, u, *args, **kwargs):

def get_fig(self): # pragma: no cover
"""
Get a figure suitable to plot the solution of this problem
Get a figure suitable to plot the solution of this problem.
Returns
-------
Expand All @@ -110,16 +121,16 @@ def get_fig(self): # pragma: no cover

def plot(self, u, t=None, fig=None, comp='u'): # pragma: no cover
r"""
Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
Plot the solution.
Parameters
----------
u : dtype_u
Solution to be plotted
t : float
Time to display at the top of the figure
fig : matplotlib.pyplot.figure.Figure
Figure with the correct structure
fig : matplotlib.pyplot.figure.Figure, optional
Figure with the same structure as a figure generated by `self.get_fig`. If none is supplied, a new figure will be generated.
Returns
-------
Expand Down Expand Up @@ -157,7 +168,18 @@ class Burgers2D(GenericSpectralLinear):
dtype_f = imex_mesh

def __init__(self, nx=64, nz=64, epsilon=0.1, fux=2, fuz=1, mode='T2U', **kwargs):
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
nx (int): Spatial resolution in x direction
nz (int): Spatial resolution in z direction
epsilon (float): viscosity
fux (int): Frequency of the initial conditions in x-direction
fuz (int): Frequency of the initial conditions in z-direction
mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
"""
self._makeAttributeAndRegister('nx', 'nz', 'epsilon', 'fux', 'fuz', 'mode', localVars=locals(), readOnly=True)

bases = [
{'base': 'fft', 'N': nx},
Expand Down
110 changes: 102 additions & 8 deletions pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,27 @@


class Heat1DChebychev(GenericSpectralLinear):
"""
1D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1) using a Chebychev spectral method.
"""

dtype_u = mesh
dtype_f = mesh

def __init__(self, nvars=128, a=0, b=0, f=1, nu=1.0, mode='T2U', **kwargs):
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
nvars (int): Resolution
a (float): Left BC value
b (float): Right BC value
f (int): Frequency of the solution
nu (float): Diffusion parameter
mode ('T2T' or 'T2U'): Mode for Chebychev method.
"""
self._makeAttributeAndRegister('nvars', 'a', 'b', 'f', 'nu', 'mode', localVars=locals(), readOnly=True)

bases = [{'base': 'chebychev', 'N': nvars}]
components = ['u', 'ux']
Expand Down Expand Up @@ -58,7 +74,18 @@ def eval_f(self, u, *args, **kwargs):
f[iu] = me[iu]
return f

def u_exact(self, t, noise=0):
def u_exact(self, t, noise=0, seed=666):
"""
Get exact solution at time `t`
Args:
t (float): When you want the exact solution
noise (float): Add noise of this level
seed (int): Random seed for the noise
Returns:
Heat1DChebychev.dtype_u: Exact solution
"""
xp = self.xp
iu, iux = self.index(self.components)
u = self.spectral.u_init
Expand All @@ -76,7 +103,7 @@ def u_exact(self, t, noise=0):
if noise > 0:
assert t == 0
_noise = self.u_init
rng = self.xp.random.default_rng(seed=666)
rng = self.xp.random.default_rng(seed=seed)
_noise[iu] = rng.normal(size=u[iu].shape)
noise_hat = self.transform(_noise)
low_pass = self.get_filter_matrix(axis=0, kmax=self.nvars - 2)
Expand All @@ -95,11 +122,25 @@ def u_exact(self, t, noise=0):


class Heat1DUltraspherical(GenericSpectralLinear):
"""
1D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1) using an ultraspherical spectral method.
"""

dtype_u = mesh
dtype_f = mesh

def __init__(self, nvars=128, a=0, b=0, f=1, nu=1.0, **kwargs):
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
nvars (int): Resolution
a (float): Left BC value
b (float): Right BC value
f (int): Frequency of the solution
nu (float): Diffusion parameter
"""
self._makeAttributeAndRegister('nvars', 'a', 'b', 'f', 'nu', localVars=locals(), readOnly=True)

bases = [{'base': 'ultraspherical', 'N': nvars}]
components = ['u']
Expand Down Expand Up @@ -147,7 +188,18 @@ def eval_f(self, u, *args, **kwargs):
f[iu][...] = me[iu]
return f

def u_exact(self, t, noise=0):
def u_exact(self, t, noise=0, seed=666):
"""
Get exact solution at time `t`
Args:
t (float): When you want the exact solution
noise (float): Add noise of this level
seed (int): Random seed for the noise
Returns:
Heat1DUltraspherical.dtype_u: Exact solution
"""
xp = self.xp
iu = self.index('u')
u = self.spectral.u_init
Expand All @@ -161,7 +213,7 @@ def u_exact(self, t, noise=0):
if noise > 0:
assert t == 0
_noise = self.u_init
rng = self.xp.random.default_rng(seed=666)
rng = self.xp.random.default_rng(seed=seed)
_noise[iu] = rng.normal(size=u[iu].shape)
noise_hat = self.transform(_noise)
low_pass = self.get_filter_matrix(axis=0, kmax=self.nvars - 2)
Expand All @@ -178,13 +230,34 @@ def u_exact(self, t, noise=0):


class Heat2DChebychev(GenericSpectralLinear):
"""
2D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1)x(-1,1) using spectral methods based on FFT and Chebychev.
"""

dtype_u = mesh
dtype_f = mesh

def __init__(self, nx=128, ny=128, base_x='fft', base_y='chebychev', a=0, b=0, c=0, fx=1, fy=1, nu=1.0, **kwargs):
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
nx (int): Resolution in x-direction
ny (int): Resolution in y-direction
base_x (str): Spectral base in x-direction
base_y (str): Spectral base in y-direction
a (float): BC at y=0 and x=-1
b (float): BC at y=0 and x=+1
c (float): BC at y=1 and x = -1
fx (int): Horizontal frequency of initial conditions
fy (int): Vertical frequency of initial conditions
nu (float): Diffusion parameter
"""
assert nx % 2 == 1 or base_x == 'fft'
assert ny % 2 == 1 or base_y == 'fft'
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
self._makeAttributeAndRegister(
'nx', 'ny', 'base_x', 'base_y', 'a', 'b', 'c', 'fx', 'fy', 'nu', localVars=locals(), readOnly=True
)

bases = [{'base': base_x, 'N': nx}, {'base': base_y, 'N': ny}]
components = ['u', 'ux', 'uy']
Expand Down Expand Up @@ -265,13 +338,34 @@ def u_exact(self, t):


class Heat2DUltraspherical(GenericSpectralLinear):
"""
2D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1)x(-1,1) using spectral methods based on FFT and Gegenbauer.
"""

dtype_u = mesh
dtype_f = mesh

def __init__(
self, nx=128, ny=128, base_x='fft', base_y='ultraspherical', a=0, b=0, c=0, fx=1, fy=1, nu=1.0, **kwargs
):
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
nx (int): Resolution in x-direction
ny (int): Resolution in y-direction
base_x (str): Spectral base in x-direction
base_y (str): Spectral base in y-direction
a (float): BC at y=0 and x=-1
b (float): BC at y=0 and x=+1
c (float): BC at y=1 and x = -1
fx (int): Horizontal frequency of initial conditions
fy (int): Vertical frequency of initial conditions
nu (float): Diffusion parameter
"""
self._makeAttributeAndRegister(
'nx', 'ny', 'base_x', 'base_y', 'a', 'b', 'c', 'fx', 'fy', 'nu', localVars=locals(), readOnly=True
)

bases = [{'base': base_x, 'N': nx}, {'base': base_y, 'N': ny}]
components = ['u']
Expand Down
47 changes: 43 additions & 4 deletions pySDC/implementations/problem_classes/RayleighBenard.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,30 @@

class RayleighBenard(GenericSpectralLinear):
"""
Rayleigh-Benard Convection is a variation of incompressible Navier-Stokes. See, for instance https://doi.org/10.1007/s00791-020-00332-3.
Rayleigh-Benard Convection is a variation of incompressible Navier-Stokes.
The equations we solve are
u_x + v_z = 0
T_t - kappa (T_xx + T_zz) = -uT_x - vT_z
u_t - nu (u_xx + u_zz) + p_x = -uu_x - vu_z
v_t - nu (v_xx + v_zz) + p_z - T = -uv_x - vv_z
with u the horizontal velocity, v the vertical velocity (in z-direction), T the temperature, p the pressure, indices
denoting derivatives, kappa=(Rayleigh * Prandl)**(-1/2) and nu = (Rayleigh / Prandl)**(-1/2). Everything on the left
hand side, that is the viscous part, the pressure gradient and the buoyancy due to temperature are treated
implicitly, while the non-linear convection part on the right hand side is integrated explicitly.
The domain, vertical boundary conditions and pressure gauge are
Omega = [0, 8) x (-1, 1)
T(z=+1) = 0
T(z=-1) = 2
u(z=+-1) = v(z=+-1) = 0
integral over p = 0
The spectral discretization uses FFT horizontally, implying periodic BCs, and an ultraspherical method vertically to
facilitate the Dirichlet BCs.
Parameters:
Prandl (float): Prandl number
Expand All @@ -29,13 +52,25 @@ def __init__(
self,
Prandl=1,
Rayleigh=2e6,
nx=257,
nx=256,
nz=64,
BCs=None,
dealiasing=3 / 2,
comm=None,
**kwargs,
):
"""
Constructor. `kwargs` are forwarded to parent class constructor.
Args:
Prandl (float): Prandtl number
Rayleigh (float): Rayleigh number
nx (int): Resolution in x-direction
nz (int): Resolution in z direction
BCs (dict): Vertical boundary conditions
dealiasing (float): Dealiasing for evaluating the non-linear part in real space
comm (mpi4py.Intracomm): Space communicator
"""
BCs = {} if BCs is None else BCs
BCs = {
'T_top': 0,
Expand Down Expand Up @@ -108,8 +143,10 @@ def __init__(
M_lhs = {i: {i: U02 @ Id} for i in ['u', 'v', 'T']}
self.setup_M(M_lhs)

# Prepare going from second (first for divergence free equation) derivative basis back to Chebychov-T
self.base_change = self._setup_operator({**{comp: {comp: S2} for comp in ['u', 'v', 'T']}, 'p': {'p': S1}})

# BCs
self.add_BC(
component='p', equation='p', axis=1, v=self.BCs['p_integral'], kind='integral', line=-1, scalar=True
)
Expand Down Expand Up @@ -238,7 +275,7 @@ def get_fig(self): # pragma: no cover

def plot(self, u, t=None, fig=None, quantity='T'): # pragma: no cover
r"""
Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
Plot the solution.
Parameters
----------
Expand All @@ -247,7 +284,9 @@ def plot(self, u, t=None, fig=None, quantity='T'): # pragma: no cover
t : float
Time to display at the top of the figure
fig : matplotlib.pyplot.figure.Figure
Figure with the correct structure
Figure with the same structure as a figure generated by `self.get_fig`. If none is supplied, a new figure will be generated.
quantity : (str)
quantity you want to plot
Returns
-------
Expand Down
Loading

0 comments on commit 0e8e997

Please sign in to comment.