Navier-Stokes equations in cylindrical geometry

[Yelyzaveta Velizhanina]

Hello all,

I would like to set up a problem that requires solving Navier-Stokes equations in cylindrical geometry. As Dedalus does not currently provide the interface for cylindrical coordinate system, I found a possible workaround in a folder d3_scratch -- test_disk_pipe.py where you combine polar bases and 1D Fourier basis. It would be exactly what I need if it worked. When I'm trying to run this script or set up my own script in a similar fashion, I get an error:

Traceback (most recent call last):
  File "/Users/yelyzavetavelizhanina/work/subcritical/scripts/tests/pipe.py", line 89, in <module>
    solver = solvers.InitialValueSolver(problem, timesteppers.RK111)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/solvers.py", line 494, in __init__
    self.subproblems = subsystems.build_subproblems(self, self.subsystems, build_matrices=['M', 'L'])
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/subsystems.py", line 68, in build_subproblems
    build_subproblem_matrices(solver, subproblems, build_matrices)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/subsystems.py", line 80, in build_subproblem_matrices
    subproblem.build_matrices(matrices)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/subsystems.py", line 427, in build_matrices
    eqn_blocks = eqn[name].expression_matrices(subproblem=self, vars=vars, ncc_cutoff=solver.ncc_cutoff, max_ncc_terms=solver.max_ncc_terms)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/operators.py", line 711, in expression_matrices
    operator_mat = self.subproblem_matrix(subproblem)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/basis.py", line 4691, in subproblem_matrix
    matrix = U.T.conj() @ matrix
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/scipy/sparse/_base.py", line 636, in __rmatmul__
    return self._rmul_dispatch(other)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/scipy/sparse/_base.py", line 614, in _rmul_dispatch
    ret = self.transpose()._mul_dispatch(tr)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/scipy/sparse/_base.py", line 577, in _mul_dispatch
    raise ValueError('dimension mismatch')
ValueError: dimension mismatch

I noticed though that if I get rid of vector equations, the problem is built successfully. The following script:

polar = de.PolarCoordinates("phi", "r")

axial = de.Coordinate("z")

dist = de.Distributor((axial, polar), dtype=np.float64)

disk = de.DiskBasis(polar, (12, 20), dtype=np.float64, radius=1.0)

axis = de.RealFourier(axial, 10, bounds=(0, 2))

phi, r = disk.local_grids()

z = axis.local_grids()[0]

u = dist.Field(bases=(axis, disk), name="u")

v = dist.Field(bases=(axis, disk), name="v")

w = dist.Field(bases=(axis, disk), name="w")

tau_u = dist.Field(bases=(axis, disk.edge), name="tau_u")

tau_v = dist.Field(bases=(axis, disk.edge), name="tau_u")

tau_w = dist.Field(bases=(axis, disk.edge), name="tau_w")

# Substitutions

dz = lambda a: de.Differentiate(a, axial)

lap = lambda a: de.Laplacian(a, polar)

lift = lambda a: de.Lift(a, disk, -1)

# Problem

problem = de.IVP([u, v, w, tau_u, tau_v, tau_w], namespace=locals())

problem.add_equation("dt(u) - lap(u) - dz(dz(u)) + lift(tau_u) = 0")

problem.add_equation("dt(v) - lap(v) - dz(dz(v)) + lift(tau_v) = 0")

problem.add_equation("dt(w) - lap(w) - dz(dz(w)) + lift(tau_w) = -1")

problem.add_equation("u(r=1) = 0")

problem.add_equation("v(r=1) = 0")

problem.add_equation("w(r=1) = 0")

solver = problem.build_solver(de.SBDF2)

solver.stop_sim_time = 5

flow = de.GlobalFlowProperty(solver, cadence=100)

flow.add_property(u * u, name="u2")

while solver.proceed:

solver.step(1e-2)

if (solver.iteration - 1) % 100 == 0:

log_prop = np.sqrt(flow.max("u2"))

print(

f"Iteration={solver.iteration}, Time={solver.sim_time:e}, "

f"flow max={log_prop}"

)

solver.log_stats()

produces a good solution -- steady pipe flow. Once I try to set a problem where the pressure is an unknown by components, I get an error

Traceback (most recent call last):
  File "/Users/yelyzavetavelizhanina/work/subcritical/scripts/tests/test_cylinder.py", line 64, in <module>
    solver.step(1e-2)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/solvers.py", line 645, in step
    self.timestepper.step(dt, wall_elapsed)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/core/timesteppers.py", line 172, in step
    sp.LHS_solver = solver.matsolver(sp.LHS, solver)
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/dedalus/libraries/matsolvers.py", line 129, in __init__
    self.LU = spla.splu(matrix.T.tocsc(), permc_spec='NATURAL')
  File "/Users/yelyzavetavelizhanina/miniforge3/envs/dedalus/lib/python3.10/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py", line 366, in splu
    return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
RuntimeError: Factor is exactly singular

I was wondering if it's possible to set up such problem using either of these two methods?

Thanks in advance!

Best regards,

Yelyzaveta.

[Keaton Burns]

Hi Yelyzaveta,

Unfortunately this is not yet working and supported. It’s functionality that many of us are obviously very interested in, but there are some technical implementation issues that need to be addressed to get it all working.  But stay tuned — we will definitely announce this is in a new v3 release when it is added, hopefully in the next few months.

Best,

-Keaton

--
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/e8b4a2b5-2596-4275-a657-e330cd03c426n%40googlegroups.com.

[Yelyzaveta Velizhanina]

Hi Keaton,

Got it, thanks for your answer and your work :)

Best,

Yelyzaveta.

On 20 Sep 2022, at 15:02, Keaton Burns <keaton...@gmail.com> wrote:



You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/pvWx_PbrvyQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBgzKeGmuOrWBCs8Y12Z6rYVz%2BgM%3DTsYRcCb0Q6wVZ2cmQ%40mail.gmail.com.

[Ankit Barik]

Hi Keaton and other devs,

Are there any updates on this?

Cheers,

Ankit

[Keaton Burns]

Hi Ankit,

Yes we have made a lot of progress implementing cylinders this winter and have some large-scale examples running. There is still some work left to get everything working with general equations, and we expect to have that done and released later this spring.

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/68e6839d-0d1f-4ad2-8be8-803f8fdd327bn%40googlegroups.com.