Ill-posed problem

[Pierre-Yves Goffin]

Hi everyone,

I'm getting the following error when trying to run my simulation: "ValueError: Pencil (np.int64(0),) has 5 constant equations for 0 constant variables plus 3 differential equations / tau terms.", indicating that I'm not considering the right number of equations, that I use wrong boundary conditions or that my problem is ill-posed in general. I read several posts in this group, but couldn't find any way to get my things working.

I'm solving a 2D Newtonian incompressible channel flow. My variables are thus the pressure (p), the velocity components along the ex and ey directions (u,v) and a constant forcing term to impose a constant mass flow rate (Fdp). My equations are the equations for the first order formulation (I'm using Dedalus v2), the continuity equation, the momentum equations and the conditions for imposing a constant Fdp and a constant mass flow rate. The boundary conditions are a no-slip condition for the velocity and a gauge condition for the pressure.

The twist is that I'm not considering a straight channel. In fact, the width of my channel varies with x and I thus tried to use a coordinate mapping between the real coordinates (x,y) and the mapped coordinates (\xi,\eta) \in [-1,1] x [-1,1]. The first variable is just a scaled version of x, but \eta is some kind of scaled y, with the varying channel half width: \eta = y/\delta.

I would really appreciate some help on that subject. Thank you in advance for the support you provide and your time.

Pierre-Yves Goffin

P.S.: I provide an example of the problem definition below. To be clear, delta_1 and delta_2 are respectively the 1st and 2nd derivatives of \delta along \xi.

import dedalus.public as de

from dedalus.extras import flow_tools

import numpy as np

xi_basis = de.Fourier('xi',64,interval=(-1.0,1.0),dealias=3/2)

eta_basis = de.Chebyshev('eta',64,interval=(-1.0,1.0),dealias=3/2)

domain = de.Domain([xi_basis,eta_basis],grid_dtype=np.float64)

problem = de.IVP(domain,variables=['Fdp','p','u','v','u_eta','v_eta'])

problem.parameters['Re'] = Re

problem.parameters['Lx'] = Lx

problem.parameters['pi'] = np.pi

# Definition of delta

delta_max = 1.0

delta_min = 0.9

c0 = delta_min/delta_max

c1 = (delta_max-delta_min)/delta_max

problem.parameters['c0'] = c0

problem.parameters['c1'] = c1

problem.substitutions['delta'] = 'c1*(sin(pi*xi/2))**2 + c0'

problem.substitutions['delta_1'] = 'c1*sin(pi*xi)'

problem.substitutions['delta_2'] = '2*c1*cos(pi*xi)'

problem.substitutions['Adv(A,A_eta)'] = 'u*dxi(A)/Lx - delta_1*eta*u*A_eta/delta + v*A_eta/delta'

problem.substitutions['LapLinear(A)'] = 'dxi(dxi(A))/Lx**2'

problem.substitutions['LapNonLinear(A,A_eta)'] = 'deta(A_eta)/delta**2 + ((delta_1*eta/delta)**2)*deta(A_eta) - 2*eta*delta_1*dxi(A_eta)/(delta*Lx) + eta*(2*delta_1**2 - delta_2*delta)*A_eta/(delta**2)'

# 1ST ORDER FORMULATION

problem.add_equation("u_eta - deta(u) = 0")

problem.add_equation("v_eta - deta(v) = 0")

# CONTINUITY

problem.add_equation("dxi(u)/Lx = -v_eta/delta + delta_1*eta*u_eta/delta")

# MOMENTUM

problem.add_equation("dt(u) + dxi(p)/Lx - LapLinear(u)/Re + Fdp = delta_1*eta*deta(p)/delta - Adv(u,u_eta) + LapNonLinear(u,u_eta)/Re")

problem.add_equation("dt(v) - LapLinear(v)/Re = -deta(p)/delta - Adv(v,v_eta) + LapNonLinear(v,v_eta)/Re")

# FORCING

problem.add_equation("Fdp = 0", condition="(nxi != 0)")

problem.add_equation("deta(Fdp) = 0", condition="(nxi == 0)")

problem.add_equation("integ(u,'eta') = 2/delta", condition="(nxi == 0)")

# BOUNDARY CONDITIONS

problem.add_bc('left(u) = 0')

problem.add_bc('right(u) = 0')

problem.add_bc('left(v) = 0')

problem.add_bc('right(v) = 0', condition="(nxi != 0)")

problem.add_bc('left(p) = 0', condition="(nxi == 0)")

# BUILD THE SOLVER

solver = problem.build_solver(eval('de.timesteppers.' + 'RK443'))

[Keaton Burns]

Hi Pierre-Yves,

Dedalus isn’t recognizing the LHS of the momentum equations are differential equations in eta because they don’t contain any eta derivatives. The best approach is to probably linearize the RHS Laplacian term and place the leading order constant coefficient term (corresponding to the mean channel height, or something) on the LHS, so that the equations have the same linear structure as for a normal channel, and just with corrections on the RHS.

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 visit https://groups.google.com/d/msgid/dedalus-users/0fb29ab9-6019-4fdd-9cc5-6f9c69295b35n%40googlegroups.com.

[Pierre-Yves Goffin]

Hi Keaton,

First of all, thank you for the quick answer.

I'm think I understand what you mean when you say that the LHS is missing some eta derivatives, but I'm not sure to understand your point when it gets to linearizing the Laplacian of the RHS. Do you mean that I should linearize the non-constant coefficients of the RHS differential operators ?

I tried that with the first non-constant coefficient (first term of "LapNonLinear"). I replaced deta(A_eta)/delta**2 with (3-2*delta)*deta(A_eta), where the function 1/delta**2 is replaced by its linear approximation around delta ≃ 1. Then I added the 3*deta(A_eta) term in the "LapLinear" substitution and modified the first term of the "LapNonLinear" substitution to -2*delta*deta(A_eta):

problem.substitutions['LapLinear(A,A_eta)'] = 'dxi(dxi(A))/L**2 + 3*deta(A_eta)'

problem.substitutions['LapNonLinear(A,A_eta)'] = '-2*delta*deta(A_eta) + ((delta_1*eta/delta)**2)*deta(A_eta) - 2*eta*delta_1*dxi(A_eta)/(delta*L) + eta*(2*delta_1**2 - delta_2*delta)*A_eta/(delta**2)'

However I must admit that I'm a bit puzzled with that method as the Laplacian is thus not what it should be as I linearized one of the coefficients. Am I right ?

Anyway, I tried that and now I'm getting a "Factor is exactly singular" error. So I imagine that I understood you wrong.

Pierre-Yves

[Keaton Burns]

Hi Pierre-Yves,

I haven’t gone through your particular notation in detail, but the idea is just to have the average part (if the channel didn’t vary) of the operators on the LHS, and only the deviations on the RHS. So if you put the channel deviations to zero in your equations, the extra terms in the RHS would be zero and the system would be identical to the normal channel case.

Best,

-Keaton

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/c55acae5-d4db-4e08-8b18-abc0255e9a5dn%40googlegroups.com.

[Pierre-Yves Goffin]

Hi Keaton,

Ok thank you I understand better what you meant.

Here I thus simply decomposed the term deta(A_eta)/delta**2 into deta(A_eta) + (1/delta**2 - 1)*deta(A_eta). The second term can thus go to the RHS (and vanishes when delta=1) and the first term is the classical wall-normal second derivative, which can go to the LHS.

problem.substitutions['LapLinear(A,A_eta)'] = 'dxi(dxi(A))/L**2 + deta(A_eta)'

problem.substitutions['LapNonLinear(A,A_eta)'] = '(1/delta**2 - 1)*deta(A_eta) + ((delta_1*eta/delta)**2)*deta(A_eta) - 2*eta*delta_1*dxi(A_eta)/(delta*L) + eta*(2*delta_1**2 - delta_2*delta)*A_eta/(delta**2)'

However I'm still getting the same "Factor is exactly singular" error, so I don't know if there is another issue somewhere else of if that did not solve everything.

Pierre-Yves

[Keaton Burns]

I think you need to linearize/separate all of the terms with delta in this fashion, so that the leading order parts are on the LHS everywhere as they are in the normal channel problem, for instance in the divergence equation as well.

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/c6b03fd8-1846-4793-990a-751cbf315375n%40googlegroups.com.

[Pierre-Yves Goffin]

Hi Keaton,

Thank you very much, it worked !

In the meantime I had tried doing that with the continuity equation too, but that didn't help. What I was missing was the pressure gradient term in the momentum equation, whose LHS didn't have the same shape as the one of the straight channel. I thus applied that for the pressure gradient term too and it worked.

Thank you for your support and your reactivity.

Regards,
Pierre-Yves