Questions about tau method and difficulty creating square system
Ford Khoudary
import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as d3
from scipy.special import jv
import logging
logger = logging.getLogger('__name__')
#Parameters
Nphi = 32
Nr = 128
l = 0.7
D = 1*10^(-7)
a = 100
b = 1*10^(-5)
s = 5*10^(-6)
stop_sim_time = 20
timestep = 1*10^(-3)
timestepper = d3.RK222
dtype = np.float64
#Bases
coords = d3.PolarCoordinates('phi','r')
dist = d3.Distributor(coords, dtype=dtype)
disk = d3.DiskBasis(coords, shape=(Nphi, Nr), radius=1, dtype=dtype)
edge = disk.edge
#Fields
p = dist.Field(name='p', bases=disk) #density
u = dist.VectorField(coords, name='u', bases=disk) #velocity
tau_u = dist.VectorField(coords, name='tau_u', bases=disk.edge)
tau_p = dist.Field(name='tau_p ')
t = dist.Field()
#Substitutions
phi, r = dist.local_grids(disk)
#Initial Conditions
p0 = 0.1
p['g'][0]=p0
u.fill_random('g', seed=42, distribution='normal')
#Problem
problem = d3.IVP([p, u], time=t, namespace=locals())
problem.add_equation("dt(p)=-div(p*u)")
problem.add_equation("dt(u)+s*grad(p)-D*lap(u)= (-l*u)@grad(u)+(a-b*(u@u))*u")
problem.add_equation("radial(u(r=1))=0")
#Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time
# Main loop
try:
logger.info('Starting main loop')
while solver.proceed:
solver.step(timestep)
if (solver.iteration-1) % 100 == 0:
max_u = 'u'
logger.info("Iteration=%i, Time=%e, dt=%e, max(u)=%e" %(solver.iteration, solver.sim_time, timestep, max_u))
except:
logger.error('Exception raised, triggering end of main loop.')
raise
finally:
solver.log_stats()
Jeffrey S. Oishi
I'm glad you're excited about using Dedalus! Your script has a number
of issues, some of which come from a lack of familiarity with python
syntax:
D = 1*10^(-7)
is not what you want. You want
D = 1e-7.
Carat (^) does not do exponentiation in python; I suggest you read a
basic python tutorial to get your head around syntactic stuff like
this. This is a dangerous problem because carat is indeed a
well-defined python operator, so it does not raise a syntax error.
Instead it silently gives you a very, very different numerical value.
There are also some mathematical issues with your formulation; I have
fixed those here: with diffusive hydrodynamics (one laplacian) in the
disk, you get one boundary condition and thus one (vector) tau because
there is one laplacian and one boundary. That tau goes on the left
hand side of the u equation. Since the density field has no
laplacians, it has no boundary conditions.
The attached script runs, though I think the domain might be too small
to see much action. I've also attached a very simple plotting script.
I've added some output that might be useful to you, including the
magnitude of p = u/umag as defined in the paper and a few other
things. I've also corrected the fact that in the paper you pointed out
alpha and beta are density-dependent, rendering the entire generalized
drag term non-linear. If alpha is a constant, then that part of the
drag becomes linear and can be moved to the LHS.
The most important thing is boundary conditions: right now, it is easy
to do no-slip (u(r=1) = 0) in the disk but not so easy to do other
conditions. If you need different boundary conditions, let us know; I
know you specified no normal flow. You actually need more than that
for viscous flow; you would need to specify the other component of
velocity for the system to be well-posed.
Finally, I'm very interested in generalized hydrodynamics and active
matter (and have worked a bit with Toner-Tu before), so I'm very
excited to see what you find with this! I'm also happy to chat
off-line if you like.
Jeff
> You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/8aac6b24-2474-4694-8ef3-bd86fcb8bf4fn%40googlegroups.com.