Implement @ operator for dot product

Author: kburns

README.md

@@ -27,7 +27,7 @@ For example, to simulate incompressible hydrodynamics in a ball, you can symboli
 
 ```python
 problem.add_equation("div(u) + tau_p = 0")
-problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - dot(u,grad(u))")
+problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - u@grad(u)")
 problem.add_equation("u(r=1) = 0")
 problem.add_equation("integ(p) = 0")
 ```

dedalus/core/field.py

@@ -110,6 +110,16 @@ def __rmul__(self, other):
         from .arithmetic import Multiply
         return Multiply(other, self)
 
+    def __matmul__(self, other):
+        # Call: self @ other
+        from .arithmetic import DotProduct
+        return DotProduct(self, other)
+
+    def __rmathmul__(self, other):
+        # Call: other @ self
+        from .arithmetic import DotProduct
+        return DotProduct(other, self)
+
     def __truediv__(self, other):
         # Call: self / other
         return (self * other**(-1))

docs/pages/changes_from_v2.rst

@@ -96,7 +96,7 @@ Vector calculus operators
 Along with vector and tensor-valued fields, vectorial differential operators (``Gradient``, ``Divergence``, ``Curl``, and ``Laplacian``) are now available.
 This dramatically simplifies the symbolic specification of vector and tensor-valued equations, particularly in curvilinear coordinates.
 Individual partial derivative operators are now usually just used in 1D equations.
-Vector algebra operations (dot products, cross products, and outer products) are also available through the ``Dot``, ``Cross``, and regular multiplication operators.
+Vector algebra operations (dot products, cross products, and outer products) are also available through the ``Dot()`` or ``@``, ``Cross()``, and ``*`` operators.
 
 For instance, an operator for computing the strain rate tensor from a velocity field can be created like:
 

docs/pages/gauge_conditions.rst

@@ -18,7 +18,7 @@ At first glance, we might think to build the problem variables and specify the e
     # Problem
     problem = d3.IVP([p, u], namespace=locals())
     problem.add_equation("div(u) = 0")
-    problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - dot(u,grad(u))")
+    problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - u@grad(u)")
 
 This formulation produces square matrices (same number of modes in the equations as the variables), but they are not all nonsingular.
 The problem is that the pressure gauge is undetermined -- any constant can be added to the pressure without affecting any of the equations.
@@ -42,7 +42,7 @@ Another option is to expand the system by adding a spatially-constant variable (
     # Problem
     problem = d3.IVP([p, u, tau_p], namespace=locals())
     problem.add_equation("div(u) + tau_p = 0")
-    problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - dot(u,grad(u))")
+    problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - u@grad(u)")
     problem.add_equation("integ(p) = 0")
 
 We've added a single additional degree of freedom to the variables and a single additional constraint, so this system is still square.

examples/ivp_2d_rayleigh_benard/rayleigh_benard.py

@@ -71,8 +71,8 @@
 # First-order form: "lap(f)" becomes "div(grad_f)"
 problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
 problem.add_equation("trace(grad_u) + tau_p = 0")
-problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - dot(u,grad(b))")
-problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - dot(u,grad(u))")
+problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")
+problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - u@grad(u)")
 problem.add_equation("b(z=0) = Lz")
 problem.add_equation("u(z=0) = 0")
 problem.add_equation("b(z=Lz) = 0")
@@ -100,7 +100,7 @@
 
 # Flow properties
 flow = d3.GlobalFlowProperty(solver, cadence=10)
-flow.add_property(np.sqrt(d3.dot(u,u))/nu, name='Re')
+flow.add_property(np.sqrt(u@u)/nu, name='Re')
 
 # Main loop
 startup_iter = 10

examples/ivp_2d_shear_flow/shear_flow.py

@@ -56,8 +56,8 @@
 
 # Problem
 problem = d3.IVP([u, s, p, tau_p], namespace=locals())
-problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - dot(u,grad(u))")
-problem.add_equation("dt(s) - D*lap(s) = - dot(u,grad(s))")
+problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - u@grad(u)")
+problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")
 problem.add_equation("div(u) + tau_p = 0")
 problem.add_equation("integ(p) = 0") # Pressure gauge
 
@@ -87,7 +87,7 @@
 
 # Flow properties
 flow = d3.GlobalFlowProperty(solver, cadence=10)
-flow.add_property(d3.dot(u,ez)**2, name='w2')
+flow.add_property((u@ez)**2, name='w2')
 
 # Main loop
 try:

examples/ivp_ball_internally_heated_convection/internally_heated_convection.py

@@ -80,7 +80,7 @@
 problem = d3.IVP([p, u, T, tau_p, tau_u, tau_T], namespace=locals())
 problem.add_equation("div(u) + tau_p = 0")
 problem.add_equation("dt(u) - nu*lap(u) + grad(p) - r_vec*T + lift(tau_u) = - cross(curl(u),u)")
-problem.add_equation("dt(T) - kappa*lap(T) + lift(tau_T) = - dot(u,grad(T)) + kappa*T_source")
+problem.add_equation("dt(T) - kappa*lap(T) + lift(tau_T) = - u@grad(T) + kappa*T_source")
 problem.add_equation("shear_stress = 0")  # Stress free
 problem.add_equation("radial(u(r=1)) = 0")  # No penetration
 problem.add_equation("radial(grad(T)(r=1)) = -2")
@@ -117,7 +117,7 @@
 
 # Flow properties
 flow = d3.GlobalFlowProperty(solver, cadence=10)
-flow.add_property(d3.dot(u,u), name='u2')
+flow.add_property(u@u, name='u2')
 
 # Main loop
 try:

examples/ivp_disk_libration/libration.py

@@ -66,7 +66,7 @@
 # Problem
 problem = d3.IVP([p, u, tau_u, tau_p], time=t, namespace=locals())
 problem.add_equation("div(u) + tau_p = 0")
-problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - dot(u, grad(u0)) - dot(u0, grad(u))")
+problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - u@grad(u0) - u0@grad(u)")
 problem.add_equation("u(r=1) = 0")
 problem.add_equation("integ(p) = 0")
 
@@ -82,11 +82,11 @@
 snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.1, max_writes=20)
 snapshots.add_task(-d3.div(d3.skew(u)), scales=(4, 1), name='vorticity')
 scalars = solver.evaluator.add_file_handler('scalars', sim_dt=0.01)
-scalars.add_task(d3.integ(0.5*d3.dot(u,u)), name='KE')
+scalars.add_task(d3.integ(0.5*u@u), name='KE')
 
 # Flow properties
 flow = d3.GlobalFlowProperty(solver, cadence=100)
-flow.add_property(d3.dot(u,u), name='u2')
+flow.add_property(u@u, name='u2')
 
 # Main loop
 try:

examples/ivp_shell_convection/shell_convection.py

@@ -73,8 +73,8 @@
 # Problem
 problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
 problem.add_equation("trace(grad_u) + tau_p = 0")
-problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - dot(u,grad(b))")
-problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*er + lift(tau_u2) = - dot(u,grad(u))")
+problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")
+problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*er + lift(tau_u2) = - u@grad(u)")
 problem.add_equation("b(r=Ri) = 1")
 problem.add_equation("u(r=Ri) = 0")
 problem.add_equation("b(r=Ro) = 0")
@@ -91,10 +91,7 @@
 b['g'] += (Ri - Ri*Ro/r) / (Ri - Ro) # Add linear background
 
 # Analysis
-db = b# - d3.Average(b, coords.S2coordsys)
-flux = d3.dot(er, -kappa*d3.grad(db) + db*u)
-#db.store_last = True
-flux.store_last = True
+flux = er @ (-kappa*d3.grad(b) + u*b)
 snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=10, max_writes=10)
 snapshots.add_task(b(r=(Ri+Ro)/2), scales=dealias, name='bmid')
 snapshots.add_task(flux(r=Ro), scales=dealias, name='flux_r_outer')
@@ -109,7 +106,7 @@
 
 # Flow properties
 flow = d3.GlobalFlowProperty(solver, cadence=10)
-flow.add_property(np.sqrt(d3.dot(u,u))/nu, name='Re')
+flow.add_property(np.sqrt(u@u)/nu, name='Re')
 
 # Main loop
 try:

examples/ivp_sphere_shallow_water/shallow_water.py

@@ -65,7 +65,7 @@
 # Initial conditions: balanced height
 c = dist.Field(name='c')
 problem = d3.LBVP([h, c], namespace=locals())
-problem.add_equation("g*lap(h) + c = - div(dot(u, grad(u)) + 2*Omega*zcross(u))")
+problem.add_equation("g*lap(h) + c = - div(u@grad(u) + 2*Omega*zcross(u))")
 problem.add_equation("ave(h) = 0")
 solver = problem.build_solver()
 solver.solve()
@@ -79,7 +79,7 @@
 
 # Problem
 problem = d3.IVP([u, h], namespace=locals())
-problem.add_equation("dt(u) + nu*lap(lap(u)) + g*grad(h) + 2*Omega*zcross(u) = - dot(u, grad(u))")
+problem.add_equation("dt(u) + nu*lap(lap(u)) + g*grad(h) + 2*Omega*zcross(u) = - u@grad(u)")
 problem.add_equation("dt(h) + nu*lap(lap(h)) + H*div(u) = - div(h*u)")
 
 # Solver