Skip to content

Commit

Permalink
Adopts new approach of ncc_cutoff as solver rather than problem keyword.
Browse files Browse the repository at this point in the history
  • Loading branch information
bpbrown committed Apr 15, 2021
1 parent f8b983b commit e21d186
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 6 deletions.
8 changes: 4 additions & 4 deletions mdwarf_hydro.py
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ def source_function(r):
Phi = trace(dot(e, e)) - 1/3*(trace_e*trace_e)

#Problem
problem = problems.IVP([u, p, s, τ_u, τ_s], ncc_cutoff=ncc_cutoff)
problem = problems.IVP([u, p, s, τ_u, τ_s])
problem.add_equation((ddt(u) + grad(p) - Co2*T*grad(s) - Ek*ρ_inv*viscous_terms + LiftTau(τ_u,-1),
- dot(u, e) - cross(ez_g, u)), condition = "ntheta != 0")
#problem.add_equation((ρ*ddt(u) + ρ*grad(p) - Co2*ρ*T*grad(s) - Ek*viscous_terms + LiftTau(τ_u,-1),
Expand All @@ -243,11 +243,11 @@ def source_function(r):

if args['--thermal_equilibrium']:
logger.info("solving for thermal equilbrium")
equilibrium = problems.LBVP([s, τ_s], ncc_cutoff=ncc_cutoff)
equilibrium = problems.LBVP([s, τ_s])
equilibrium.add_equation((-(lap(s)+ dot(grad_lnT, grad(s))) + LiftTau(τ_s,-1),
ρ*source))
equilibrium.add_equation((s(r=radius), 0))
eq_solver = solvers.LinearBoundaryValueSolver(equilibrium)
eq_solver = solvers.LinearBoundaryValueSolver(equilibrium, ncc_cutoff=ncc_cutoff)
eq_solver.solve()

amp = 1e-2
Expand All @@ -269,7 +269,7 @@ def source_function(r):
s['g'] += amp*noise

# Solver
solver = solvers.InitialValueSolver(problem, timesteppers.SBDF2)
solver = solvers.InitialValueSolver(problem, timesteppers.SBDF2, ncc_cutoff=ncc_cutoff)
timestepper_history = [0,1]

reducer = GlobalArrayReducer(d.comm_cart)
Expand Down
4 changes: 2 additions & 2 deletions structure.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,13 @@ def lane_emden(Nmax, Lmax=0, m=1.5, n_rho=3, radius=1,
lap = lambda A: operators.Laplacian(A, c)
Pow = lambda A,n: operators.Power(A,n)
LiftTau = lambda A: operators.LiftTau(A, br, -1)
problem = problems.NLBVP([f, R, τ], ncc_cutoff=ncc_cutoff)
problem = problems.NLBVP([f, R, τ])
problem.add_equation((lap(f) + LiftTau(τ), - R**2 * Pow(f,m)))
problem.add_equation((f(r=0), 1))
problem.add_equation((f(r=radius), np.exp(-n_rho/m, dtype=dtype))) # explicit typing to match domain

# Solver
solver = solvers.NonlinearBoundaryValueSolver(problem)
solver = solvers.NonlinearBoundaryValueSolver(problem, ncc_cutoff=ncc_cutoff)
# Initial guess
f['g'] = np.cos(np.pi/2 * r)**2
R['g'] = 5
Expand Down

0 comments on commit e21d186

Please sign in to comment.