Problem 3D Elasticity Equation; Non-square system:

[DedELas]

Dear Dedalus-users,

I am trying to implement a PDE in three dimensions. Following the old standard of D2, first I broke down the PDE to a first order-system and introduced appropriate tau-terms. However, when I try to solve it, the error "ValueError: Non-square system: group=(0, 1, None), I=156, J=144" occurred. In 2D it's working fine, however in 3D something goes wrong. I could see any error. Thank you for your help!

Enclosed you can find the code for two (working) and three dimension (not working). I am using Dedalus v3, it's installed by conda.

import numpy as np
import dedalus.public as d3

# Parameters
Lx, Ly = 2*np.pi, np.pi
Nx, Ny = 12, 12
dtype = np.float64

E = 100
nu = 0.3
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))

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

# Fields
u1 = dist.Field(name='u1', bases=(xbasis, ybasis))
u2 = dist.Field(name='u2', bases=(xbasis, ybasis))
u1x = dist.Field(name='u1x', bases=(xbasis, ybasis))
u1y = dist.Field(name='u1y', bases=(xbasis, ybasis))
u2x = dist.Field(name='u2x', bases=(xbasis, ybasis))
u2y = dist.Field(name='u2y', bases=(xbasis, ybasis))
tau_1 = dist.Field(name='tau_1', bases=xbasis)
tau_2 = dist.Field(name='tau_2', bases=xbasis)
tau_3 = dist.Field(name='tau_3', bases=xbasis)
tau_4 = dist.Field(name='tau_4', bases=xbasis)

# Substitutions
dy = lambda A: d3.Differentiate(A, coords['y'])
dx = lambda A: d3.Differentiate(A, coords['x'])
lift_basis = ybasis.derivative_basis(1)
lift = lambda A, n: d3.Lift(A, lift_basis, n)

fla = dist.Field()
fmu = dist.Field()
fla['g'] = lmbda
fmu['g'] = mu

problem = d3.LBVP([u1, u2, u1y, u2y, tau_1, tau_2, tau_3, tau_4], namespace=locals())

problem.add_equation("dx((2*fmu+fla)*dx(u1)) + dy(fla*u1y) + dx(fla*u2y) + dy(fla*dx(u2)) + lift(tau_1,-1)= 0")
problem.add_equation("dy(fla*dx(u1)) + dx(fla*u1y) + dy((2*fmu+fla)*u2y) + dx(fla*dx(u2)) + lift(tau_3,-1)= 0")
problem.add_equation("u1y - dy(u1) + lift(tau_2,-1) = 0")
problem.add_equation("u2y - dy(u2) + lift(tau_4,-1) = 0")

problem.add_equation("u1(y=0) = 0")
problem.add_equation("u1(y=Ly) = 4")
problem.add_equation("u2(y=0) = 0")
problem.add_equation("u2(y=Ly) = 2")

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

# -----------------------------------------------

#

# -----------------------------------------------

3D Code

import numpy as np
import dedalus.public as d3

# Parameters
Lx, Ly, Lz = 2, 2, 2
Nx, Ny, Nz = 12, 12, 12
dtype = np.float64

E = 100
nu = 0.3
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))

# Bases
coords = d3.CartesianCoordinates('x', 'y', 'z')
dist   = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(0, Ly))
zbasis = d3.Chebyshev(coords['z'], size=Nz, bounds=(0, Lz))

# Fields
u1 = dist.Field(name='u1', bases=(xbasis, ybasis, zbasis))
u2 = dist.Field(name='u2', bases=(xbasis, ybasis, zbasis))
u3 = dist.Field(name='u3', bases=(xbasis, ybasis, zbasis))

u1z = dist.Field(name='u1z', bases=(xbasis, ybasis, zbasis))
u2z = dist.Field(name='u2z', bases=(xbasis, ybasis, zbasis))
u3z = dist.Field(name='u3z', bases=(xbasis, ybasis, zbasis))

tau01 = dist.Field(name='tau01', bases=xbasis)
tau02 = dist.Field(name='tau02', bases=xbasis)
tau03 = dist.Field(name='tau03', bases=xbasis)
tau04 = dist.Field(name='tau04', bases=xbasis)
tau05 = dist.Field(name='tau05', bases=xbasis)
tau06 = dist.Field(name='tau06', bases=xbasis)


# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
dz = lambda A: d3.Differentiate(A, coords['z'])
lift_basis_z = zbasis.derivative_basis(1)
lift_z = lambda A, n: d3.Lift(A, lift_basis_z, n)

fla = dist.Field()
fmu = dist.Field()
fla['g'] = 0.
fmu['g'] = mu

problem = d3.LBVP([u1,u2,u3, u1z,u2z,u3z, tau01, tau02, tau03, tau04, tau05, tau06], namespace=locals())

problem.add_equation("dx((2*fmu+fla)*dx(u1)) + dy(fmu*(dy(u1) + dx(u2))) + dz(fmu*(u1z + dx(u3))) \
+ lift_z(tau01,-1) = 0")

problem.add_equation("dx(fmu*(dy(u1) + dx(u2))) + dz(fmu*(u2z + dy(u3))) + dy((2*fmu+fla)*dy(u2))\
+ lift_z(tau02,-1) = 0")

problem.add_equation("dx(fmu*(u1z + dx(u3))) + dy(fmu*(u2z + dy(u3))) + dz((2*fmu+fla)*u3z)\
+ lift_z(tau03,-1) = 0")

problem.add_equation("u1z - dz(u1) + lift_z(tau04,-1) = 0")
problem.add_equation("u2z - dz(u2) + lift_z(tau05,-1) = 0")
problem.add_equation("u3z - dz(u3) + lift_z(tau06,-1) = 0")

problem.add_equation("u1(z=0)  = 0")
problem.add_equation("u1(z=Ly) = 4")
problem.add_equation("u2(z=0)  = 0")
problem.add_equation("u2(z=Lz) = 0")
problem.add_equation("u3(z=0)  = 0")
problem.add_equation("u3(z=Lz) = 0")

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

[Daniel Lecoanet]

Hi,

For the 3D problem, I think your tau variables need to be defined over the x and y bases. These terms are essentially tau(x, y)*T(z), where T is a Chebyshev polynomial.

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/e93478c5-37ac-431f-ab5d-e9303b79aa7bn%40googlegroups.com.

[DedELas]

Hi Daniel,

thank you for help. You are right, it's working, if use Chebyshev for z-basis and Fourier for y- and x-basis. However, if I use Chebyshev polynominal for y- and z-basis and Fourier for x-basis,  I got strange and wrong results. It's also unclear to me, how to select degree of Chebyshev Polynomial for the tau terms?

Please find below the adapted code with Fourier discretization in x-direction and Chebyshev in y- and z-direction.

---------------------------------------------------------

import numpy as np
import dedalus.public as d3

from pyevtk.hl import gridToVTK

# Parameters
Lx, Ly, Lz = 2, 2, 2

Nx, Ny, Nz = 24, 12, 12

dtype = np.float64

E = 100
nu = 0.3
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))

# Bases
coords = d3.CartesianCoordinates('x','y','z')
dist   = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))

ybasis = d3.Chebyshev(coords['y'], size=Ny, bounds=(0, Ly))

zbasis = d3.Chebyshev(coords['z'], size=Nz, bounds=(0, Lz))

# Fields
u1 = dist.Field(name='u1', bases=(xbasis, ybasis, zbasis))
u2 = dist.Field(name='u2', bases=(xbasis, ybasis, zbasis))
u3 = dist.Field(name='u3', bases=(xbasis, ybasis, zbasis))

u1y = dist.Field(name='u1y', bases=(xbasis, ybasis, zbasis))
u2y = dist.Field(name='u2y', bases=(xbasis, ybasis, zbasis))
u3y = dist.Field(name='u3y', bases=(xbasis, ybasis, zbasis))

u1z = dist.Field(name='u1z', bases=(xbasis, ybasis, zbasis))
u2z = dist.Field(name='u2z', bases=(xbasis, ybasis, zbasis))
u3z = dist.Field(name='u3z', bases=(xbasis, ybasis, zbasis))

tau01 = dist.Field(name='tau01', bases=(xbasis,ybasis))
tau02 = dist.Field(name='tau02', bases=(xbasis,ybasis))
tau03 = dist.Field(name='tau03', bases=(xbasis,ybasis))
tau04 = dist.Field(name='tau04', bases=(xbasis,ybasis))
tau05 = dist.Field(name='tau05', bases=(xbasis,ybasis))
tau06 = dist.Field(name='tau06', bases=(xbasis,ybasis))
tau07 = dist.Field(name='tau07', bases=(xbasis,ybasis))
tau08 = dist.Field(name='tau08', bases=(xbasis,ybasis))
tau09 = dist.Field(name='tau09', bases=(xbasis,ybasis))
tau10 = dist.Field(name='tau10', bases=(xbasis,ybasis))
tau11 = dist.Field(name='tau11', bases=(xbasis,ybasis))
tau12 = dist.Field(name='tau12', bases=(xbasis,ybasis))

# Substitutions
dx = lambda A: d3.Differentiate(A, coords['x'])
dy = lambda A: d3.Differentiate(A, coords['y'])
dz = lambda A: d3.Differentiate(A, coords['z'])

lift_basis_y = ybasis.derivative_basis(1)
lift_y = lambda A, n: d3.Lift(A, lift_basis_z, n)

lift_basis_z = zbasis.derivative_basis(1)
lift_z = lambda A, n: d3.Lift(A, lift_basis_z, n)

fla = dist.Field()
fmu = dist.Field()

fla['g'] = lmbda
fmu['g'] = mu

problem = d3.LBVP([u1,u2,u3, u1y,u2y,u3y, u1z,u2z,u3z, tau01, tau02, tau03, tau04, tau05, tau06, tau07\
, tau08, tau09, tau10, tau11, tau12], namespace=locals())

problem.add_equation("dx((2*fmu+fla)*dx(u1) + dx(fla*u2y) + dx(fla*u3z)) + dy(fmu*(u1y + dx(u2))) + \
dz(fmu*(u1z + dx(u3))) + lift_z(tau01,-1) + lift_y(tau02,-3) = 0")

problem.add_equation("dx(fmu*(u1y + dx(u2))) + dy(fla*dx(u1) + dy((2*fmu+fla)*u2y) + dy(fla*u3z)) +\
dz(fmu*(u2z + u3y)) + lift_z(tau03,-1) + lift_y(tau04,-3) = 0")

problem.add_equation("dx(fmu*(u1z + dx(u3))) + dy(fmu*(u2z + u3y)) + \
dz(fla*dx(u1) + fla*u2y + (2*fmu+fla)*u3z) + lift_z(tau05,-1) + lift_y(tau06,-3) = 0")

problem.add_equation("u1z - dz(u1) + lift_z(tau07,-5) = 0")
problem.add_equation("u2z - dz(u2) + lift_z(tau08,-5) = 0")
problem.add_equation("u3z - dz(u3) + lift_z(tau09,-5) = 0")
problem.add_equation("u1y - dy(u1) + lift_z(tau10,-5) = 0")
problem.add_equation("u2y - dy(u2) + lift_z(tau11,-5) = 0")
problem.add_equation("u3y - dy(u3) + lift_z(tau12,-5) = 0")

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

problem.add_equation("u1(y=Ly) = 0")

problem.add_equation("u1(z=0)  = 0")

problem.add_equation("u1(z=Lz) = 0")

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

problem.add_equation("u2(y=Ly) = 0")

problem.add_equation("u2(z=0)  = 0")

problem.add_equation("u2(z=Lz) = 1")

problem.add_equation("u3(y=0)  = 0")
problem.add_equation("u3(y=Ly) = 0")

problem.add_equation("u3(z=0)  = 0")
problem.add_equation("u3(z=Lz) = 0")


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

u1g = u1.allgather_data('g')
u2g = u2.allgather_data('g')
u3g = u3.allgather_data('g')

# ---------------- export results to paraview -------------------------
# Dimensions
dx, dy, dz = Lx/Nx, Ly/Ny, Lz/Nz
ncells = Nx * Ny * Nz
# Coordinates
X = np.arange(0, Lx + 0.1*dx, dx, dtype='float64')
Y = np.arange(0, Ly + 0.1*dy, dy, dtype='float64')
Z = np.arange(0, Lz + 0.1*dz, dz, dtype='float64')
x = np.zeros((Nx+1, Ny+1, Nz+1))
y = np.zeros((Nx+1, Ny+1, Nz+1))
z = np.zeros((Nx+1, Ny+1, Nz+1))
#
for k in range(Nz+1):
    for j in range(Ny + 1):
        for i in range(Nx + 1):
            x[i,j,k] = X[i] + dx
            y[i,j,k] = Y[j] + dy
            z[i,j,k] = Z[k] + dz
           
           
gridToVTK("./Results", x, y, z, cellData = {
    "u1" : u1g.reshape((Nx,Ny,Nz)),
    "u2" : u2g.reshape((Nx,Ny,Nz)),
    "u3" : u3g.reshape((Nx,Ny,Nz))})

[Daniel Lecoanet]

Hi,

I don’t think multiple Chebyshev bases are implemented right now. Also, the tau terms for them are complicated due to corner conditions.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/03a0af26-6ec0-4b38-99db-68cf04719732n%40googlegroups.com.