Adding info about Robin and similar boundary conditions

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)

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