ValueError: Problem is coupled along distributed dimensions:

[rumayel pallock]

Hello Dedalus Users,

I am trying to solve a 3D LBVP problem where the domain is periodic in x, rho and dirichlet in y.

My goal is to run this in MPI. I have attached two files. When i try to run LBVP_Not_Working.py using 4 cores, i am getting following error:

"ValueError: Problem is coupled along distributed dimensions: (0, 1)"

These are the equations I am using,

# Linear operator equation

LU_l = (1/Re)*(-dx(p1) + 2*dx(dx(U1)) + dy(U1y) + dx(V1y) - RaEr*(dx(Q11_1)) - RaEr*(dy(Q12_1))) + lift(tau_U12)

LU_r =  - (U1*dx(U) + V1*dy(U)) - (U*dx(U1) + V*dy(U1))

LV_l = (1/Re)*(-dy(p1) + dx(U1y) + dx(dx(V1)) + 2*dy(V1y) - RaEr*(dx(Q12_1) - dy(Q11_1))) + lift(tau_V12)

LV_r =  - (U1*dx(V) + V1*dy(V)) - (U*dx(V1) + V*dy(V1))

LQ11_l = lam*dx(U1) + Q11_1  + dx(dx(Q11_1)) + dy(Q11_1y)  + lift(tau_Q11_21)

LQ11_r = - (30*Q11*Q11*Q11_1 + 10*Q12*Q12*Q11_1 + 20*Q12*Q12_1*Q11)- (U*dx(Q11_1) + V*dy(Q11_1) + U1*dx(Q11) + V1*dy(Q11)) + (dy(U))*Q12_1 + (dy(U1))*Q12 - (dx(V))*Q12_1 - (dx(V1))*Q12

LQ12_l = (lam/2)*(dy(U1) + dx(V1)) + Q12_1  + dx(dx(Q12_1)) + dy(Q12_1y)  + lift(tau_Q12_21)

LQ12_r = - (10*Q11**2*Q12_1 + 30*Q12**2*Q12_1 + 20*Q11*Q12*Q11_1)- (U*dx(Q12_1) + V*dy(Q12_1) + U1*dx(Q12) + V1*dy(Q12)) - (dy(U))*Q11_1 - (dy(U1))*Q11 + (dx(V))*Q11_1 + (dx(V1))*Q11

# Divergence free condition

divU1 = (dx(U1)) + (V1y) + tau_p1

# Define N: Residuals

N1 = (1/T)*drho(U) - FU

N2 = (1/T)*drho(V) - FV

N3 = (1/T)*drho(Q11) - FQ11

N4 = (1/T)*drho(Q12) - FQ12

G2_rhs = integral((1/T**2)*(drho(U)*N1 + drho(V)*N2 + drho(Q11)*N3 + drho(Q12)*N4))

Equation_LHS_U = (1/T)*drho(U1) - LU_l - LU_r

Equation_RHS_U = -(1/T)*drho(U) + FU + (1/T**2)*G2*drho(U)

Equation_LHS_V = (1/T)*drho(V1) - LV_l - LV_r

Equation_RHS_V = -(1/T)*drho(V) + FV + (1/T**2)*G2*drho(V)

Equation_LHS_Q11 = (1/T)*drho(Q11_1) - LQ11_l - LQ11_r

Equation_RHS_Q11 = -(1/T)*drho(Q11) + FQ11 + (1/T**2)*G2*drho(Q11)

Equation_LHS_Q12 = (1/T)*drho(Q12_1) - LQ12_l - LQ12_r

Equation_RHS_Q12 = -(1/T)*drho(Q12) + FQ12 + (1/T**2)*G2*drho(Q12)

# problem

problem = d3.LBVP([U1,V1,Q11_1,Q12_1,p1,tau_U11,tau_U12,tau_V11,tau_V12,tau_Q11_11,tau_Q11_21,tau_Q12_11,tau_Q12_21,tau_p1], namespace = locals())

problem.add_equation("Equation_LHS_U = Equation_RHS_U")

problem.add_equation("Equation_LHS_V = Equation_RHS_V")

problem.add_equation("Equation_LHS_Q11 = Equation_RHS_Q11")

problem.add_equation("Equation_LHS_Q12 = Equation_RHS_Q12")

problem.add_equation("divU1 = 0")

problem.add_equation("integ_xy(p1) = 0") # pressure gauge

problem.add_equation("U1(y=0)=0")

problem.add_equation("U1(y=h)=0")

problem.add_equation("V1(y=0)=0")

problem.add_equation("V1(y=h)=0")

problem.add_equation("Q11_1(y=0) = 0")

problem.add_equation("Q11_1(y=h)= 0")

problem.add_equation("Q12_1(y=0) = 0")

problem.add_equation("Q12_1(y=h) = 0")

Now if i delete LU_r, LV_r, LQ11_r, LQ12_r then I can run in MPI without error. LBVP_Working.py is the version where those are deleted.

Does anyone know how I can keep those terms and run the program in MPI? Any help will be appreciated.

Regards,

Rumayel

[Daniel Lecoanet]

Hi Rumayel,

To solve a linear boundary value problem, Dedalus turns your equation into the form L.X = b, where X is the state vector of all of your variables, and L.X represents all terms on the LHS of the equals sign in your equations. Currently, the L matrix includes non-constant coefficients such as U which depend on all three variables x, rho, and y. That means that the problem is fully coupled, i.e., there is a single L the problem. In contrast, if the NCCs were independent of x, then each Fourier mode in x would decouple from each other, and you would have Nx different L matrices (which depend on the wavenumber). The latter case can be parallelized: each core constructs and performs matrix solves for some subset of L matrices. However, since your problem is fully-coupled, we do not allow the matrix construction/matrix solve of the single L matrix to be parallelized. If you wanted to do the matrix solve in parallel, you could get L and b from Dedalus, and then pass it to your favorite parallel matrix solve library.

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 visit https://groups.google.com/d/msgid/dedalus-users/f9b7ea78-8e23-4f1d-bc62-4fc620c07177n%40googlegroups.com.
<LBVP_Not_Working.py><LBVP_Working.py>

[rumayel pallock]

Hi Daniel,

Thanks for your reply.

Let me know if I have misunderstood or not. Because my equations have NCC(U,V,Q11,Q12...), so that I cannot construct the matrix using multiple core because they are coupled. Instead I have to use a single core to construct the matrix and then I can use other libraries to solve in parallel.

I have a following question. For my problem, if Nx= 16, Ny = 8 and Nrho = 16, dedalus can build the matrix and solve in single core.

But If i increase the number of modes, (e.g: changing Nrho = 16 to 32), i get "killed" message. My RAM is 32GB. Even if i try to run in computers with Higher memory in cluster, i get the same message. Do you know or can you share any resources on how to tackle this kind of situation? 

Rumayel

[rumayel pallock]

Also can you tell me how can I access L and b from dedalus? 

Regards,

Rumayel

[Daniel Lecoanet]

I believe your problem, being fully-coupled, has only a single subproblem. Then I think the L matrix is given by

solver.subproblems[0].L_min

And the b vector is constructed via

for field in solver.F:

    field.change_layout('c')

solver.subproblems[0].gather_outputs(solver.F)

You might be interested in looking at the code for how we solve LBVPs here:

Daniel

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/6266e973-5965-48ee-9608-dcc411b231c7n%40googlegroups.com.

[rumayel pallock]

Hi Daniel,

Thanks for your reply. I really appreciate your help.

I am getting this error when i tried to get L_min: AttributeError: 'Subproblem' object has no attribute 'L_min'

So I have checked all the attributes by using following code at the end:

# Build solver

timestep = 1e-6

solver = problem.build_solver()

attributes = dir(solver.subproblems[0])

print(attributes)

and these are the available attributes:

['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_build_buffer_views', '_compressed_buffer', '_input_buffer', '_input_field_views', '_input_views', '_output_buffer', '_output_field_views', '_output_views', 'build_matrices', 'check_condition', 'coeff_shape', 'coeff_size', 'coeff_slices', 'dist', 'domain', 'dtype', 'expand_matrices', 'field_shape', 'field_size', 'field_slices', 'gather_inputs', 'gather_outputs', 'group', 'group_dict', 'inclusion_matrices', 'problem', 'scatter_inputs', 'scatter_outputs', 'shape', 'solver', 'subsystems', 'valid_modes']

There is not attribute named "L_min" for LBVP solver.

In your reference link, I saw L_min inside "EigenvalueSolver". But did not find anything like that in LBVP. Let me know If I am missing something.

Rumayel

[Daniel Lecoanet]

Oh, you might have to build the L matrix first. Try:

solver.build_matrices(solver.subproblems[0], ['L'])

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/ea4f7708-80bd-4159-86eb-09886f531335n%40googlegroups.com.

[rumayel pallock]

Hi Daniel,

I could have build the L matrix and RHS. 

Now I have a different problem. I have moved to a simple poison problem with a NCC. The equation is 

NCC*laplacian(u(x,y)) + tau_p = f

and, integ(u) = 0; 

where f is the forcing term, and the problem is periodic in both x and y.

I have kept the grid size Nx = 32 and Ny = 32.

When I check the size of the matrix, i expect the matrix to be 1025x1025 (32*32+1 = 1025)

But actually, the matrix size is 962x962 (31*31+1 = 962)

Now my question is, why is it N-1 for fourier basis? is it the row corresponding to -sin(0x) is removed?

I have added the code for reference.

[Daniel Lecoanet]

Yes, that’s right, the -sin(0x) component is dropped. This is done by marking that mode as invalid (https://github.com/DedalusProject/dedalus/blob/c829c76bb2d2a896ebd7905fab4a7c7be60e941d/dedalus/core/basis.py#L1123) which is then used in the preconditioners when constructing the matrices.

Daniel

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/4d5e6e23-0e70-4af4-9d60-1d8f2645def8n%40googlegroups.com.
<sample.py>

[rumayel pallock]

Hi Daniel,

I have a question regarding RHS of the equations.

For the following case:

Tx = dx(T) + lift(tau_T1)

equation_rhs = dx(f) + f

# define problem

problem = d3.LBVP([T,tau_T1,tau_T2],namespace = locals())

problem.add_equation("dx(Tx) + lift(tau_T2) = equation_rhs")

problem.add_equation("T(x=0) = 0")

problem.add_equation("T(x=L) = 0")

if i do,

solver = problem.build_solver()

# Find RHS using built-in method

# LHS matrix

solver.build_matrices(solver.subproblems, ['L'])

csr = solver.subproblems[0].L_min

# Compute RHS

solver.evaluator.evaluate_scheduled(iteration = solver.iteration)

for field in solver.F:

    field.change_layout('c')

RHS = solver.subproblems[0].gather_outputs(solver.F)

print(RHS)

It gives me the RHS. But I should also recreate the RHS from following way,

equation_rhs_ev = equation_rhs.evaluate()

RHS_manual = np.concatenate([equation_rhs_ev['c'].ravel(), np.zeros(2)])  # add two zeros for the two BCs

print(RHS_manual)

For some reason, RHS and RHS_manual are not the same if i stick to 'c' (I mean field.change_layout('c') and equation_rhs_ev['c']).  

Otherwise, if I use 'g' for both of them, they give me same vector. Can you tell me why for 'c', I am not getting same? Is there something that I need to do? I am trying to create the RHS without building the solver as for big problem, it is taking time.

I have attached the full code for your ease.

Regards,

Rumayel

[Daniel Lecoanet]

Hi Rumayel,

They are different because you are calculating the coefficients of different bases. Your by-hand calculation is in the Jacobi basis with a=b=1/2, whereas the Dedalus calculation is using the Jacobi basis with a=b=5/2. You should probably be using a=b=3/2, which you would get if you set the lift_basis to derivative_basis(1). In any case, when you convert to the same basis, you appear to get the same numbers. See attached script.

Daniel

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/32ee50e7-fdf1-4eaa-8c88-8af43488675en%40googlegroups.com.
<small_problem.py>