search

DedalusProject / dedalus

A flexible framework for solving PDEs with modern spectral methods.
http://dedalus-project.org/
GNU General Public License v3.0
489 stars 115 forks source link

Overflow in 2D, flat and periodic active nematic simulation #263

Closed PeterHampshire closed 1 year ago

PeterHampshire commented 1 year ago

Hi thanks for creating such a nice resource. I would like to reproduce the active nematic simulations from the following paper https://arxiv.org/abs/2306.15352. The simulations are on a 2D, flat and periodic domain and the dimensionless equations are as follows: image

I attempted to solve the equations using Dedalus 3 but I get the error "RuntimeWarning: overflow encountered in multiply" on the first time step (see first code below).

Perhaps the problem arises because of the constraint equation (Equation 2 in the image above). It is linear in the velocity field but the velocity terms have coefficients that are dependent on the other variables (the nematic q tensor and density rho). There are many ways to write it such that the LHS is linear in the problem variables, so my attempt was to simply divide by rho.

My question is whether it is possible to solve these equations with Dedalus? If so, how can I change my code so that I do not get an overflow error on the first time step?

I tried reducing the max_time_step and the standard deviation of the noise in the initial conditions for q and rho but in both cases I still got an overflow error on the first time step. I also tried solving a LBVP to get a consistent initialisation for the velocity field (see second code below) but I get the error "f"Problem is coupled along distributed dimensions: {tuple(np.where(coupled_nonlocal)[0])}"".

Thanks, Peter


import numpy as np
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)

# Parameters
Nx, Ny = 256, 256
dealias = 3/2
timestepper = d3.RK222
max_timestep = 1e-4
Nsteps = 10000
dN = 100                
k_d=0.1
L=1
b=20
eta_rot=1
beta=-0.2
lambdabar_sun=6
lambdabar=1.3*(1+2*np.sqrt(k_d))**2
kappa=-0.8
c=4
a=1/2*(lambdabar_sun+c)
Lx, Ly = 8*2*np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2)), 8*2*np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2))

stop_sim_time=Nsteps*max_timestep
sim_dt=dN*max_timestep

dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=dealias)
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)

# Fields
rho = dist.Field(name='rho', bases=(xbasis,ybasis))
v = dist.VectorField(coords, name='v', bases=(xbasis,ybasis))
q = dist.TensorField((coords, coords), name='q', bases=(xbasis,ybasis))

# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
x, y = dist.local_grids(xbasis, ybasis)
ex, ey = coords.unit_vector_fields(dist)
identity = ex*ex+ey*ey
d = 1/2*(d3.grad(v) + d3.trans(d3.grad(v)))
w = 1/2*(d3.grad(v) - d3.trans(d3.grad(v)))
#qhat=d3.dt(q)+v@d3.grad(q)-w@q+q@w
d_dev = d-d3.trace(d)/2*identity
sigma = rho*(2*(d+d3.trace(d)*identity)\
              +beta*1/eta_rot*(-beta*d_dev-(2*a+b*2*d3.trace(q@q))*q+L*(d3.lap(q)+1/rho*d3.grad(q)@d3.grad(rho))+lambdabar_sun*rho*q)\
                  +lambdabar*(identity+kappa*q)\
                      -L*(2*(d3.grad(q@ex@ex)*d3.grad(q@ex@ex)+d3.grad(q@ex@ey)*d3.grad(q@ex@ey))\
                          -1/rho*q@d3.div(rho*d3.grad(q))+1/rho*d3.trans(q@d3.div(rho*d3.grad(q)))))

# Problem
problem = d3.IVP([rho, v, q], namespace=locals())
problem.add_equation("dt(rho)-lap(rho)+k_d*rho = k_d - div(rho*v)")
problem.add_equation("v = 1/rho*div(sigma)")
problem.add_equation("eta_rot*dt(q)+beta*d_dev+2*a*q-L*lap(q)=\
                     -eta_rot*(v@grad(q)-w@q+q@w)\
                         -b*2*trace(q@q)*q+L/rho*grad(q)@grad(rho)+lambdabar_sun*rho*q")

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time

# Initial conditions
rho.fill_random('g', seed=42, distribution='normal', scale = 0.05)
rho += 1
q.fill_random('g', seed=42, distribution='normal', scale = 0.05)
q['g'][1,0]=q['g'][0,1]
q['g'][1,1]=-q['g'][0,0]

# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=sim_dt, max_writes=10)
snapshots.add_task(rho, name='rho')
snapshots.add_task(d3.div(v), name='divergence')
snapshots.add_task(-d3.div(d3.skew(v)), name='vorticity')

# CFL
CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=0.2, threshold=0.1, max_change=1.5, min_change=0.5, max_dt=max_timestep)

CFL.add_velocity(v)

# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=10)

# Main loop
try:
    logger.info('Starting main loop')
    while solver.proceed:
        timestep = CFL.compute_timestep()
        solver.step(timestep)
        if (solver.iteration-1) % 100 == 0:
            #max_w = np.sqrt(flow.max('w2'))
            logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))
except:
    logger.error('Exception raised, triggering end of main loop.')
    raise
finally:
    solver.log_stats()


import numpy as np
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)

# Parameters
Nx, Ny = 256, 256
dealias = 3/2
Nsteps = 10000
dN = 100                
k_d=0.1
L=1
b=20
eta_rot=1
beta=-0.2
lambdabar_sun=6
lambdabar=1.3*(1+2*np.sqrt(k_d))**2
kappa=-0.8
c=4
a=1/2*(lambdabar_sun+c)
Lx, Ly = 8*2*np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2)), 8*2*np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2))

dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=dealias)
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)

# Fields
rho = dist.Field(name='rho', bases=(xbasis,ybasis))
v = dist.VectorField(coords, name='v', bases=(xbasis,ybasis))
q = dist.TensorField((coords, coords), name='q', bases=(xbasis,ybasis))

#Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
x, y = dist.local_grids(xbasis, ybasis)
ex, ey = coords.unit_vector_fields(dist)
identity = ex*ex+ey*ey
d = 1/2*(d3.grad(v) + d3.trans(d3.grad(v)))
w = 1/2*(d3.grad(v) - d3.trans(d3.grad(v)))
#qhat=d3.dt(q)+v@d3.grad(q)-w@q+q@w
d_dev = d-d3.trace(d)/2*identity
sigma = rho*(2*(d+d3.trace(d)*identity)\
              +beta*1/eta_rot*(-beta*d_dev-(2*a+b*2*d3.trace(q@q))*q+L*(d3.lap(q)+1/rho*d3.grad(q)@d3.grad(rho))+lambdabar_sun*rho*q)\
                  +lambdabar*(identity+kappa*q)\
                      -L*(2*(d3.grad(q@ex@ex)*d3.grad(q@ex@ex)+d3.grad(q@ex@ey)*d3.grad(q@ex@ey))\
                          -1/rho*q@d3.div(rho*d3.grad(q))+1/rho*d3.trans(q@d3.div(rho*d3.grad(q)))))

# Problem
problem = d3.LBVP([v], namespace=locals())
problem.add_equation("v-1/rho*div(rho*(2*(d+d3.trace(d)*identity)-beta*1/eta_rot*beta*d_dev))=\
                     1/rho*div(rho*(beta*1/eta_rot*(-(2*a+b*2*d3.trace(q@q))*q+L*(d3.lap(q)+1/rho*d3.grad(q)@d3.grad(rho))+lambdabar_sun*rho*q)\
                         +lambdabar*(identity+kappa*q)\
                             -L*(2*(d3.grad(q@ex@ex)*d3.grad(q@ex@ex)+d3.grad(q@ex@ey)*d3.grad(q@ex@ey))\
                                 -1/rho*q@d3.div(rho*d3.grad(q))+1/rho*d3.trans(q@d3.div(rho*d3.grad(q))))))")

# Initial conditions
rho.fill_random('g', seed=42, distribution='normal', scale = 0.05)
rho += 1
q.fill_random('g', seed=42, distribution='normal', scale = 0.05)
q['g'][1,0]=q['g'][0,1]
q['g'][1,1]=-q['g'][0,0]

#LBVP
solver = problem.build_solver()
solver.solve()
bpbrown commented 1 year ago

Peter, There's a good chance that this overflow is coming about from the 1/rho term, but maybe somewhere else. If it's the 1/rho term, I've had a lot of success reformulating problems like this by solving for ln(rho) as the variable rather than rho.

For example, grad(rho)/rho -> grad(ln_rho). And likewise for 1/rho * dt(rho) -> dt(ln_rho).

The price is that you now have nonlinear terms like exp(ln_rho) anywhere you would have had rho terms previously.

But this approach can work very well if you need to keep things like rho bounded from reaching (or crossing) non-physical zero values. --Ben

On Thu, Aug 24, 2023 at 3:22 AM PeterHampshire @.***> wrote:

Hi thanks for creating such a nice resource. I would like to reproduce the active nematic simulations from the following paper https://arxiv.org/abs/2306.15352 https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Farxiv.org%2Fabs%2F2306.15352&data=05%7C01%7CBenjamin.P.Brown%40Colorado.EDU%7C0d7d3171ff004aa8e28c08dba48c0786%7C3ded8b1b070d462982e4c0b019f46057%7C1%7C0%7C638284693565303629%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=UnN5AWFvlY7prnl4ymTRcBNsANqWSNJTz3RkuWhvbcU%3D&reserved=0. The simulations are on a 2D, flat and periodic domain and the dimensionless equations are as follows: [image: image] https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F87777339%2F262931649-bda4c8df-9c23-4c2c-96ef-09db7eae048e.png&data=05%7C01%7CBenjamin.P.Brown%40Colorado.EDU%7C0d7d3171ff004aa8e28c08dba48c0786%7C3ded8b1b070d462982e4c0b019f46057%7C1%7C0%7C638284693565303629%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=%2BNmFM3dsGDMr8f7tilYaG1kx1ThpaVWThxMlvwjF0ek%3D&reserved=0

I attempted to solve the equations using Dedalus 3 but I get the error "RuntimeWarning: overflow encountered in multiply" on the first time step (see first code below).

Perhaps the problem arises because of the constraint equation (Equation 2 in the image above). It is linear in the velocity field but the velocity terms have coefficients that are dependent on the other variables (the nematic q tensor and density rho). There are many ways to write it such that the LHS is linear in the problem variables, so my attempt was to simply divide by rho.

My question is whether it is possible to solve these equations with Dedalus? If so, how can I change my code so that I do not get an overflow error on the first time step?

I tried reducing the max_time_step and the standard deviation of the noise in the initial conditions for q and rho but in both cases I still got an overflow error on the first time step. I also tried solving a LBVP to get a consistent initialisation for the velocity field (see second code below) but I get the error "f"Problem is coupled along distributed dimensions: {tuple(np.where(coupled_nonlocal)[0])}"".

Thanks, Peter

import numpy as np import dedalus.public as d3 import logging logger = logging.getLogger(name)

Parameters

Nx, Ny = 256, 256 dealias = 3/2 timestepper = d3.RK222 max_timestep = 1e-4 Nsteps = 10000 dN = 100 k_d=0.1 L=1 b=20 eta_rot=1 beta=-0.2 lambdabar_sun=6 lambdabar=1.3(1+2np.sqrt(k_d))*2 kappa=-0.8 c=4 a=1/2(lambdabar_sun+c) Lx, Ly = 82np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2)), 82np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2))

stop_sim_time=Nstepsmax_timestep sim_dt=dNmax_timestep

dtype = np.float64

Bases

coords = d3.CartesianCoordinates('x', 'y') dist = d3.Distributor(coords, dtype=dtype) xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=dealias) ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)

Fields

rho = dist.Field(name='rho', bases=(xbasis,ybasis)) v = dist.VectorField(coords, name='v', bases=(xbasis,ybasis)) q = dist.TensorField((coords, coords), name='q', bases=(xbasis,ybasis))

Substitutions

dx = lambda A: d3.Differentiate(A, coords['x']) dy = lambda A: d3.Differentiate(A, coords['y']) x, y = dist.local_grids(xbasis, ybasis) ex, ey = coords.unit_vector_fields(dist) identity = exex+eyey d = 1/2(d3.grad(v) + d3.trans(d3.grad(v))) w = 1/2(d3.grad(v) - d3.trans(d3.grad(v))) @.**@*.**@*. d_dev = d-d3.trace(d)/2identity sigma = rho(2(d+d3.trace(d)identity)\ @*.**@*.(rho))+lambdabar_sunrhoq)\ +lambdabar(identity+kappaq)\ @.@@.@@.@@.@ey))\ @*.**@*.**(rhod3.grad(q)))))

Problem

problem = d3.IVP([rho, v, q], namespace=locals()) problem.add_equation("dt(rho)-lap(rho)+k_drho = k_d - div(rhov)") problem.add_equation("v = 1/rhodiv(sigma)") problem.add_equation("eta_rotdt(q)+betad_dev+2aq-Llap(q)=\ @.**@*.**@*.)\ @*.**@*.**(rho)+lambdabar_sunrho*q")

Solver

solver = problem.build_solver(timestepper) solver.stop_sim_time = stop_sim_time

Initial conditions

rho.fill_random('g', seed=42, distribution='normal', scale = 0.05) rho += 1 q.fill_random('g', seed=42, distribution='normal', scale = 0.05) q['g'][1,0]=q['g'][0,1] q['g'][1,1]=-q['g'][0,0]

Analysis

snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=sim_dt, max_writes=10) snapshots.add_task(rho, name='rho') snapshots.add_task(d3.div(v), name='divergence') snapshots.add_task(-d3.div(d3.skew(v)), name='vorticity')

CFL

CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=0.2, threshold=0.1, max_change=1.5, min_change=0.5, max_dt=max_timestep)

CFL.add_velocity(v)

Flow properties

flow = d3.GlobalFlowProperty(solver, cadence=10)

Main loop

try: logger.info('Starting main loop') while solver.proceed: timestep = CFL.compute_timestep() solver.step(timestep) if (solver.iteration-1) % 100 == 0:

max_w = np.sqrt(flow.max('w2'))

        logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))

except: logger.error('Exception raised, triggering end of main loop.') raise finally: solver.log_stats()

import numpy as np import dedalus.public as d3 import logging logger = logging.getLogger(name)

Parameters

Nx, Ny = 256, 256 dealias = 3/2 Nsteps = 10000 dN = 100 k_d=0.1 L=1 b=20 eta_rot=1 beta=-0.2 lambdabar_sun=6 lambdabar=1.3(1+2np.sqrt(k_d))*2 kappa=-0.8 c=4 a=1/2(lambdabar_sun+c) Lx, Ly = 82np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2)), 82np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2))

dtype = np.float64

Bases

coords = d3.CartesianCoordinates('x', 'y') dist = d3.Distributor(coords, dtype=dtype) xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=dealias) ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)

Fields

rho = dist.Field(name='rho', bases=(xbasis,ybasis)) v = dist.VectorField(coords, name='v', bases=(xbasis,ybasis)) q = dist.TensorField((coords, coords), name='q', bases=(xbasis,ybasis))

Substitutions

dx = lambda A: d3.Differentiate(A, coords['x']) dy = lambda A: d3.Differentiate(A, coords['y']) x, y = dist.local_grids(xbasis, ybasis) ex, ey = coords.unit_vector_fields(dist) identity = exex+eyey d = 1/2(d3.grad(v) + d3.trans(d3.grad(v))) w = 1/2(d3.grad(v) - d3.trans(d3.grad(v))) @.**@*.**@*. d_dev = d-d3.trace(d)/2identity sigma = rho(2(d+d3.trace(d)identity)\ @*.**@*.(rho))+lambdabar_sunrhoq)\ +lambdabar(identity+kappaq)\ @.@@.@@.@@.@ey))\ @*.**@*.**(rhod3.grad(q)))))

Problem

problem = d3.LBVP([v], namespace=locals()) problem.add_equation("v-1/rhodiv(rho(2(d+d3.trace(d)identity)-beta1/eta_rotbeta*d_dev))=\ @.**@.(rho))+lambdabar_sunrhoq)\ +lambdabar(identity+kappaq)\ @.@@.@@.@@.@ey))\ @.**@.(rho*d3.grad(q))))))")

Initial conditions

rho.fill_random('g', seed=42, distribution='normal', scale = 0.05) rho += 1 q.fill_random('g', seed=42, distribution='normal', scale = 0.05) q['g'][1,0]=q['g'][0,1] q['g'][1,1]=-q['g'][0,0]

LBVP

solver = problem.build_solver() solver.solve()

— Reply to this email directly, view it on GitHub https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FDedalusProject%2Fdedalus%2Fissues%2F263&data=05%7C01%7CBenjamin.P.Brown%40Colorado.EDU%7C0d7d3171ff004aa8e28c08dba48c0786%7C3ded8b1b070d462982e4c0b019f46057%7C1%7C0%7C638284693565303629%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=LclzWcciLxqbUvWa99dkK5waGPufJVU41Eu9d33ZOhM%3D&reserved=0, or unsubscribe https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABN7BYP6NK5KO2HTNHCNVHDXW4TORANCNFSM6AAAAAA34ZILLQ&data=05%7C01%7CBenjamin.P.Brown%40Colorado.EDU%7C0d7d3171ff004aa8e28c08dba48c0786%7C3ded8b1b070d462982e4c0b019f46057%7C1%7C0%7C638284693565303629%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=N4QFbd0Jd9EE5H2GYZITglhggLEZkwlQVZQB5NHkuAM%3D&reserved=0 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

PeterHampshire commented 1 year ago

Thanks for the suggestion @bpbrown but unfortunately I still get the same error on the first time step "RuntimeWarning: overflow encountered in multiply". Here's the code I used


#Re-writing in terms of ln(rho)

import numpy as np
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)

# Parameters
Nx, Ny = 256, 256
dealias = 3/2
timestepper = d3.RK222
max_timestep = 1e-4
Nsteps = 10000
dN = 100                
k_d=0.1
L=1
b=20
eta_rot=1
beta=-0.2
lambdabar_sun=6
lambdabar=1.3*(1+2*np.sqrt(k_d))**2
kappa=-0.8
c=4
a=1/2*(lambdabar_sun+c)
Lx, Ly = 8*2*np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2)), 8*2*np.pi/(np.sqrt(np.sqrt(k_d))/np.sqrt(2))

stop_sim_time=Nsteps*max_timestep
sim_dt=dN*max_timestep

dtype = np.float64

# Bases
coords = d3.CartesianCoordinates('x', 'y')
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(-Lx/2, Lx/2), dealias=dealias)
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)

# Fields 
lnrho = dist.Field(name='lnrho', bases=(xbasis,ybasis))
v = dist.VectorField(coords, name='v', bases=(xbasis,ybasis))
q = dist.TensorField((coords, coords), name='q', bases=(xbasis,ybasis))

# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
x, y = dist.local_grids(xbasis, ybasis)
ex, ey = coords.unit_vector_fields(dist)
identity = ex*ex+ey*ey
d = 1/2*(d3.grad(v) + d3.trans(d3.grad(v)))
w = 1/2*(d3.grad(v) - d3.trans(d3.grad(v)))
#qhat=d3.dt(q)+v@d3.grad(q)-w@q+q@w
d_dev = d-d3.trace(d)/2*identity
sigma = (2*(d+d3.trace(d)*identity)\
              +beta*1/eta_rot*(-beta*d_dev-(2*a+b*2*d3.trace(q@q))*q+L*(d3.lap(q)+d3.grad(q)@d3.grad(lnrho))+lambdabar_sun*np.exp(lnrho)*q)\
                  +lambdabar*(identity+kappa*q)\
                      -L*(2*(d3.grad(q@ex@ex)*d3.grad(q@ex@ex)+d3.grad(q@ex@ey)*d3.grad(q@ex@ey))\
                          -q@d3.div(d3.grad(q))-q@d3.grad(lnrho)@d3.grad(q)+d3.trans(q@d3.div(d3.grad(q))+q@d3.grad(lnrho)@d3.grad(q))))

# Problem
problem = d3.IVP([lnrho, v, q], namespace=locals())
problem.add_equation("dt(lnrho)-lap(lnrho)+div(v) = (dx(lnrho)**2+dy(lnrho)**2)-k_d+k_d*np.exp(-lnrho) - -v@grad(lnrho)")
problem.add_equation("v = div(sigma)+grad(lnrho)@sigma")
problem.add_equation("eta_rot*dt(q)+beta*d_dev+2*a*q-L*lap(q)=\
                     -eta_rot*(v@grad(q)-w@q+q@w)\
                         -b*2*trace(q@q)*q+L*grad(q)@grad(lnrho)+lambdabar_sun*np.exp(lnrho)*q")

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time

# Initial conditions
lnrho.fill_random('g', seed=42, distribution='normal', scale = 0.05)
#rho += 1
q.fill_random('g', seed=42, distribution='normal', scale = 0.05)
q['g'][1,0]=q['g'][0,1]
q['g'][1,1]=-q['g'][0,0]

# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=sim_dt, max_writes=10)
snapshots.add_task(lnrho, name='lnrho')
snapshots.add_task(d3.div(v), name='divergence')
snapshots.add_task(-d3.div(d3.skew(v)), name='vorticity')
#snapshots.add_task(S, name='pressure')
#snapshots.add_task(-d3.div(d3.skew(u)), name='vorticity')

# CFL
CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=0.2, threshold=0.1, max_change=1.5, min_change=0.5, max_dt=max_timestep)

CFL.add_velocity(v)

# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=10)
#flow.add_property((u@ez)**2, name='w2')

# Main loop
try:
    logger.info('Starting main loop')
    while solver.proceed:
        timestep = CFL.compute_timestep()
        solver.step(timestep)
        if (solver.iteration-1) % 100 == 0:
            #max_w = np.sqrt(flow.max('w2'))
            logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))
except:
    logger.error('Exception raised, triggering end of main loop.')
    raise
finally:
    solver.log_stats()
kburns commented 1 year ago

We use the issue tracker here for bugs and feature requests for Dedalus as a whole. I think this conversation should be continued on the mailing list, though, which is a better place for questions about formulating specific equation sets.