boundary condition depends on other variables.

[Chris Li]

Hi, I'm very new to the Dedalus project, and I'm working on solving the Navier-Stokes equations with Bousinesq approximation in an annulus geometry. The velocity field expressed as vec{u}=u \vec{e_theta} + v\vec{e_r}, and the boundary conditions are 

u = v = 0 at r=Ri and r=Ro

\frac{\partial T}{\partial r} = 0 at r=Ri

T = (1-sin(\theta))/2 at r=Ro

I followed the code in post: https://groups.google.com/g/dedalus-users/c/Zen6EhH6VEs/m/SaKi3mwGAAAJ and modified the boundary conditions as

phi, r = dist.local_grids(annulus)
sin = dist.Field(name='sin', bases=annulus.outer_edge)
sin['g']=np.sin(phi)

dr=lambda A: d3.radial(d3.grad(A)(r=Ri))

problem.add_equation("u(r=Ri) = 0")
problem.add_equation("u(r=Ro) = 0")
problem.add_equation("dr(T) = 0")
problem.add_equation("T(r=Ro) = (1-sin)/2")

Then it gives the following error:

--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[42], line 6 4 while solver.proceed: 5 timestep = CFL.compute_timestep() ----> 6 solver.step(timestep) 7 if (solver.iteration-1) % 10 == 0: 8 max_Re = flow.max('Re') File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/solvers.py:654, in InitialValueSolver.step(self, dt) 652 self.warmup_time = wall_time 653 # Advance using timestepper --> 654 self.timestepper.step(dt, wall_time) 655 # Update iteration 656 self.iteration += 1 File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/timesteppers.py:596, in RungeKuttaIMEX.step(self, dt, wall_time) 594 # Compute F(n,i-1), only doing output on first evaluation 595 if i == 1: --> 596 evaluator.evaluate_scheduled(iteration=iteration, wall_time=wall_time, sim_time=solver.sim_time, timestep=dt) 597 else: 598 evaluator.evaluate_group('F') File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/evaluator.py:92, in Evaluator.evaluate_scheduled(self, **kw) 90 """Evaluate all scheduled handlers.""" 91 handlers = [h for h in self.handlers if h.check_schedule(**kw)] ---> 92 self.evaluate_handlers(handlers, **kw) File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/evaluator.py:111, in Evaluator.evaluate_handlers(self, handlers, id, **kw) 109 fields = self.get_fields(tasks) 110 self.require_coeff_space(fields) --> 111 tasks = self.attempt_tasks(tasks, id=id) 113 # Oscillate through layouts until all tasks are evaluated 114 # Limit to 10 passes to break on potential infinite loops 115 n_layouts = len(self.dist.layouts) File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/evaluator.py:201, in Evaluator.attempt_tasks(tasks, **kw) 199 unfinished = [] 200 for task in tasks: --> 201 output = task['operator'].attempt(**kw) 202 if output is None: 203 unfinished.append(task) File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/future.py:226, in Future.attempt(self, id) 224 def attempt(self, id=None): 225 """Recursively attempt to evaluate operation.""" --> 226 return self.evaluate(id=id, force=False) File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/future.py:168, in Future.evaluate(self, id, force) 166 a.change_scales(a.domain.dealias) 167 if isinstance(a, Future): --> 168 a_eval = a.evaluate(id=id, force=force) 169 # If evaluation succeeds, substitute result 170 if a_eval is not None: File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/future.py:168, in Future.evaluate(self, id, force) 166 a.change_scales(a.domain.dealias) 167 if isinstance(a, Future): --> 168 a_eval = a.evaluate(id=id, force=force) 169 # If evaluation succeeds, substitute result 170 if a_eval is not None: File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/future.py:196, in Future.evaluate(self, id, force) 193 out.preset_scales(self.domain.dealias) 195 # Perform operation --> 196 self.operate(out) 198 # Reset to free temporary field arguments 199 self.reset() File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/arithmetic.py:250, in AddFields.operate(self, out) 248 # Set output layout 249 out.preset_layout(arg0.layout) --> 250 np.add(arg0.data, arg1.data, out=out.data) ValueError: non-broadcastable output operand with shape (1,1) doesn't match the broadcast shape (16,1)

I think I must do something wrong with the sin boundary condition, but I don't know how to fix it. Thank you so much for your help.

Best regards,

Chris Li

[Daniel Lecoanet]

Hi Chris,

I believe this is a bug we recently identified… Try (sin-1)/2 instead. Alternatively define a field, call it T_out, which is equal to (1-np.sin(phi))/2 and set T(r=To)=T_out. Let me know if those work.

Daniel

--
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/f444ab99-1168-4926-b858-1a9241abac79n%40googlegroups.com.

[Chris Li]

Hi Daniel,

Thanks for your quick reply. (sin-1)/2 do works. 

Best,

Chris

[Chris Li]

Hi Daniel,

I'm trying to understand these lines of code from your previous reply in other post https://groups.google.com/g/dedalus-users/c/Zen6EhH6VEs/m/SaKi3mwGAAAJ

# Fields p = dist.Field(name='p', bases=annulus) T = dist.Field(name='T', bases=annulus) u = dist.VectorField(coords, name='u', bases=annulus) tau_p = dist.Field(name='tau_p') tau_T1 = dist.Field(name='tau_T1', bases=annulus.inner_edge) tau_T2 = dist.Field(name='tau_T2', bases=annulus.inner_edge) tau_u1 = dist.VectorField(coords, name='tau_u1', bases=annulus.inner_edge) tau_u2 = dist.VectorField(coords, name='tau_u2', bases=annulus.inner_edge)
I'm wondering why tau_T1, tau_T2, tau_u1, and tau_u2 are using the inner_edge? How do we know which base should we use for annulus case?

Best,

Chris

On Monday, August 12, 2024 at 9:00:25 a.m. UTC-7 daniel.lecoanet wrote:

[Daniel Lecoanet]

Hi Chris,

These tau variables are functions of phi only. The annulus has an inner_edge and an outer_edge, but both are equivalent for the purpose of defining these tau variables (because you should not do anything with the tau variables that requires knowing the radius at which they are defined).

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/eb371f05-5599-4f86-9b73-da3cbc47a21fn%40googlegroups.com.

[Chris Li]

Hi Daniel,

I have another question for my system, attached in the picture. I'm trying to solve it in annulus domain, note there is a term PrRaTe_y, so I think I should convert it into e_y=sin(theta)e_r + cos(theta)e_theta first. I know that e_r can be defined as 

er = dist.VectorField(coords, bases=annulus.radial_basis)
er['g'][1] = 1

but how should I define for e_theta since I don't see there is phi_basis or theta_basis under the annulus class.

Best,

Chris

[Daniel Lecoanet]

Hi Chris,

In polar coordinates, er and etheta do not depend on r or theta. So I think you can define them as vector fields without any bases and take er[‘g’][1] = 1 and etheta[‘g’][0] = 1.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/249c4fe9-4d9b-4a67-8194-97f25f71728an%40googlegroups.com.
<system.png>

[Chris Li]

Hi Daniel,

Thanks for your reply. I'm trying to follow https://groups.google.com/g/dedalus-users/c/Zen6EhH6VEs/m/SaKi3mwGAAAJ and https://github.com/DedalusProject/dedalus/blob/master/examples/ivp_annulus_centrifugal_convection/centrifugal_convection.py to solve my system attached in the previous post. I have attached my code in the thread, and counter the following error. I don't quite understand which step I started to make error. Thanks for your help and time. I appreciate it.

--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In[10], line 21 19 while solver.proceed: 20 timestep = CFL.compute_timestep() ---> 21 solver.step(timestep) 22 if (solver.iteration-1) % 10 == 0: 23 max_Re = flow.max('Re') File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/solvers.py:654, in InitialValueSolver.step(self, dt) 652 self.warmup_time = wall_time 653 # Advance using timestepper --> 654 self.timestepper.step(dt, wall_time) 655 # Update iteration 656 self.iteration += 1 File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/core/timesteppers.py:627, in RungeKuttaIMEX.step(self, dt, wall_time) 625 else: 626 sp.LHS = (sp.M_min + k_Hii*sp.L_min) # CREATES TEMPORARY --> 627 sp.LHS_solvers[i] = solver.matsolver(sp.LHS, solver) 628 # Slice out valid subdata, skipping invalid components 629 spRHS = RHS.get_subdata(sp) File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/dedalus/libraries/matsolvers.py:141, in _SuperluFactorizedBase.__init__(self, matrix, solver) 139 elif self.trans == "H": 140 matrix = matrix.H --> 141 self.LU = spla.splu(matrix.tocsc(), 142 permc_spec=self.permc_spec, 143 diag_pivot_thresh=self.diag_pivot_thresh, 144 relax=self.relax, 145 panel_size=self.panel_size, 146 options=self.options) File ~/anaconda3/envs/dedalus3/lib/python3.12/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py:438, in splu(A, permc_spec, diag_pivot_thresh, relax, panel_size, options) 435 if (_options["ColPerm"] == "NATURAL"): 436 _options["SymmetricMode"] = True --> 438 return _superlu.gstrf(N, A.nnz, A.data, indices, indptr, 439 csc_construct_func=csc_construct_func, 440 ilu=False, options=_options) RuntimeError: Factor is exactly singular

[Daniel Lecoanet]

Hi Chris,

You’re missing tau terms. I encourage you to look at this page on recommended ways on how to set up your tau terms:

https://dedalus-project.readthedocs.io/en/latest/pages/tau_method.html

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/b51a5951-47e7-4cdd-9c8b-60b1dfaf80b0n%40googlegroups.com.
<ns_annulus.py>