From e14d9d5a63f38ffa9de90dd0a7b8411e3a974cf5 Mon Sep 17 00:00:00 2001 From: Mikael Mortensen Date: Mon, 5 Aug 2024 15:11:40 +0200 Subject: [PATCH] Adding info about Robin and similar boundary conditions --- demo/Kuramato_Sivashinsky.py | 21 ++++++++--------- docs/source/gettingstarted.rst | 41 ++++++++++++++++++++++++++-------- 2 files changed, 43 insertions(+), 19 deletions(-) diff --git a/demo/Kuramato_Sivashinsky.py b/demo/Kuramato_Sivashinsky.py index 51864a26..deeedf9d 100644 --- a/demo/Kuramato_Sivashinsky.py +++ b/demo/Kuramato_Sivashinsky.py @@ -11,14 +11,15 @@ Use Fourier basis V and tensor product space VxV """ -from sympy import symbols, exp, pi +import sys +from sympy import symbols, exp, pi, Rational import numpy as np import matplotlib.pyplot as plt from shenfun import * # Use sympy to set up initial condition x, y = symbols("x,y", real=True) -ue = exp(-0.01*((x-30*pi)**2+(y-30*pi)**2)) # + exp(-0.02*((x-15*np.pi)**2+(y)**2)) +ue = exp(-((x-30*pi)**2+(y-30*pi)**2)/100) - Rational(1, 36) / pi # Size of discretization N = (128, 128) @@ -51,8 +52,6 @@ def NonlinearRHS(self, U, U_hat, dU, gradu, **params): gradu = TVp.backward(1j*K*U_hat, gradu) dU = Tp.forward(0.5*(gradu[0]*gradu[0]+gradu[1]*gradu[1]), dU) dU.mask_nyquist(mask) - if comm.Get_rank() == 0: - dU[0, 0] = 0 return -dU #initialize @@ -61,17 +60,18 @@ def NonlinearRHS(self, U, U_hat, dU, gradu, **params): U_hat.mask_nyquist(mask) # Integrate using an exponential time integrator -def update(self, u, u_hat, t, tstep, plot_step, **params): +def update(self, u, u_hat, t, tstep, plot_step, wash, **params): if not hasattr(self, 'fig'): self.fig = plt.figure() self.cm = plt.get_cmap('hot') self.image = plt.contourf(X[0], X[1], U, 256, cmap=self.cm) - if comm.Get_rank() == 0: - u_hat[0, 0] = 0 - - if tstep % plot_step == 0 and plot_step > 0: + u_hat.mask_nyquist(mask) + if tstep % wash == 0 and wash > 0: u = u_hat.backward(u) u_hat = u.forward(u_hat) + + if tstep % plot_step == 0 and plot_step > 0: + u = u_hat.backward(u) self.image.axes.contourf(X[0], X[1], u, 256, cmap=self.cm) plt.autoscale() plt.pause(1e-6) @@ -82,10 +82,11 @@ def update(self, u, u_hat, t, tstep, plot_step, **params): if __name__ == '__main__': par = {'plot_step': 100, + 'wash': 1, 'gradu': gradu, 'count': 0} dt = 0.01 - end_time = 50 + end_time = 100 integrator = ETDRK4(T, L=LinearRHS, N=NonlinearRHS, update=update, **par) #integrator = RK4(T, L=LinearRHS, N=NonlinearRHS, update=update, **par) integrator.setup(dt) diff --git a/docs/source/gettingstarted.rst b/docs/source/gettingstarted.rst index b56864fe..dd105b7a 100644 --- a/docs/source/gettingstarted.rst +++ b/docs/source/gettingstarted.rst @@ -302,16 +302,40 @@ returns a space with 4 boundary conditions (biharmonic), where ``a`` and ``b`` are the Dirichlet and Neumann values on the left boundary, whereas ``c`` and ``d`` are the values on right. +Any kind of boundary condition may be specified. For higher order +derivatives, use the form ``bc = f"u''(-1)={a}"``, or ``bc = {'left': {'N2': a}}``, +and similar for higher order. + +It is also possible to apply Robin conditions + +.. math:: + + u(x) + c u'(x) = d + +Applying this on both sides of the domain is achieved through:: + + bc = {'left': {'R': (c, d)}, 'right': {'R': (c, d)}} + +Finally, one can use the boundary condition + +.. math:: + + u'(x) + c u''(x) = d + +Applying this on both sides of the domain:: + + bc = {'left': {'W': (c, d)}, 'right': {'W': (c, d)}} + +Currently there are no hardcoded spaces implemented for these latter two boundary +conditions, so the composition of these bases is computed using a generic approach +that is slower than for the more common Dirichlet and Neumann spaces. + The Laguerre basis is used to solve problems on the half-line :math:`x \in [0, \infty)`. For this family you can only specify boundary conditions at the left boundary. However, the Poisson equation requires only one condition, and the biharmonic problem two. The solution is automatically set to zero at :math:`x \rightarrow \infty`. -Any kind of boundary condition may be specified. For higher order -derivatives, use the form ``bc = f"u''(-1)={a}"``, or ``bc = {'left': {'N2': a}}``, -and similar for higher order. - Multidimensional problems ------------------------- @@ -496,14 +520,13 @@ Curvilinear coordinates ----------------------- Shenfun can be used to solve equations using curvilinear coordinates, like polar, cylindrical -and spherical coordinates. The feature was added April 2020, and is still rather -experimental. The curvilinear coordinates are defined by the user, who +and spherical coordinates. The curvilinear coordinates are defined by the user, who needs to provide a map, i.e., the position vector, between new coordinates and the Cartesian coordinates. The basis functions of the new coordinates need not be orthogonal, but non-orthogonal is not widely tested so use with care. -In shenfun we use non-normalized natural (covariant) basis vectors. For this -reason the equations may look a little bit different than usual. For example, -in cylindrical coordinates we have the position vector +In shenfun we use non-normalized natural (covariant) basis vectors by default. +For this reason the equations may look a little bit different than usual. For +example, in cylindrical coordinates we have the position vector .. math:: :label: eq:cylpositionvector