BC in cylinder wake forward problem

[jacktang]

Hello,

I self-train deepxde by re-producing some example codes. I've tried time-independent cylinder wake forward problem, and it works fine. Geometry and BCs are defined below:

# Geometry
farfield = dde.geometry.Rectangle([xmin, ymin], [xmax, ymax])
cylinder = dde.geometry.Disk([0.2, 0.2], 0.05)
geom     = dde.geometry.CSGDifference(farfield, cylinder)

inner_rec  = dde.geometry.Rectangle([0.1, 0.1], [0.2, 0.3])
outer_dom  = dde.geometry.CSGDifference(farfield, inner_rec)
outer_dom  = dde.geometry.CSGDifference(outer_dom, cylinder)
inner_dom  = dde.geometry.CSGDifference(inner_rec, cylinder)

# BCs
...
def boundary_cylinder(x, on_boundary):
    return on_boundary and cylinder.on_boundary(x)
...

And I also want to try time-dependent cylinder wake forward problem, and I add time domain like below

# Geometry
time = dde.geometry.TimeDomain(0, 20)
farfield = dde.geometry.Rectangle([xmin, ymin], [xmax, ymax])
cylinder = dde.geometry.Disk([0.2, 0.2], 0.05)
geom     = dde.geometry.CSGDifference(farfield, cylinder)

geom_time = dde.geometry.GeometryXTime(geom, time)

inner_rec  = dde.geometry.Rectangle([0.1, 0.1], [0.2, 0.3])
outer_dom  = dde.geometry.CSGDifference(farfield, inner_rec)
outer_dom  = dde.geometry.CSGDifference(outer_dom, cylinder)
inner_dom  = dde.geometry.CSGDifference(inner_rec, cylinder)

farfield_time = dde.geometry.GeometryXTime(farfield, time)
cylinder_time = dde.geometry.GeometryXTime(cylinder, time)
inner_dom_time = dde.geometry.GeometryXTime(inner_dom, time)
outer_dom_time = dde.geometry.GeometryXTime(outer_dom, time)

# BCs
...
def boundary_cylinder(x, on_boundary):
    return on_boundary and cylinder_time.on_boundary(x)  # <- ERROR

And I get the error: too many indices for array: array is 1-dimensional, but 2 were indexed
How can I fix this?

[lululxvi]

You may first figure out which part of the code makes the code.

[jacktang]

@lululxvi I upload the python script at the end( rename to .py).

我的主要问题是，时间无关的(稳态)NS到时间相关的NS需要做哪些调整？。稳态的NS我已经可以拿到运行结果
目前我做的调整有
1)：方程的调整，加入时间t
2)：domain的调整
3)：boundary condition的调整

但是最后报错：too many indices for array: array is 1-dimensional, but 2 were indexed。
我不清楚是我代码的问题，还是deepxde的bug

#!/usr/bin/env python
# coding: utf-8

# # Time dependent cylinder wake

import deepxde as dde
import numpy as np
from matplotlib import pyplot as plt



dde.config.set_random_seed(48)
dde.config.set_default_float('float64')

xmin, xmax = 0.0, 1.0
ymin, ymax = 0.0, 0.4

rho  = 1.0
mu   = 0.02
umax = 1.0



def navier_stokes(x, y):
    """
    System of PDEs to be minimized: incompressible Navier-Stokes equation in the
    continuum-mechanics based formulations.
    
    """
    psi, p, sigma11, sigma22, sigma12 = y[:, 0:1], y[:, 1:2], y[:, 2:3], y[:, 3:4], y[:, 4:5]
    
    u =   dde.grad.jacobian(y, x, i = 0, j = 1)
    v = - dde.grad.jacobian(y, x, i = 0, j = 0)
    
    u_x = dde.grad.jacobian(u, x, i = 0, j = 0)
    u_y = dde.grad.jacobian(u, x, i = 0, j = 1)
    u_t = dde.grad.jacobian(u, x, i = 0, j = 2)
    
    v_x = dde.grad.jacobian(v, x, i = 0, j = 0)
    v_y = dde.grad.jacobian(v, x, i = 0, j = 1)
    v_t = dde.grad.jacobian(v, x, i = 0, j = 2)
    
    sigma11_x = dde.grad.jacobian(y, x, i = 2, j = 0)
    sigma12_x = dde.grad.jacobian(y, x, i = 4, j = 0)
    sigma12_y = dde.grad.jacobian(y, x, i = 4, j = 1)
    sigma22_y = dde.grad.jacobian(y, x, i = 3, j = 1)
    
    continuumx = rho * u_t + rho * (u * u_x + v * u_y) - sigma11_x - sigma12_y
    continuumy = rho * v_t + rho * (u * v_x + v * v_y) - sigma12_x - sigma22_y
    
    constitutive1 = - p + 2 * mu * u_x - sigma11
    constitutive2 = - p + 2 * mu * v_y - sigma22
    constitutive3 = mu * (u_y + v_x) - sigma12
    constitutive4 = p + (sigma11 + sigma22) / 2
    
    return continuumx, continuumy, constitutive1, constitutive2, constitutive3, constitutive4


# Geometry defintion
time = dde.geometry.TimeDomain(0, 20)
farfield = dde.geometry.Rectangle([xmin, ymin], [xmax, ymax])
cylinder = dde.geometry.Disk([0.2, 0.2], 0.05)
geom     = dde.geometry.CSGDifference(farfield, cylinder)

geom_time = dde.geometry.GeometryXTime(geom, time)

inner_rec  = dde.geometry.Rectangle([0.1, 0.1], [0.2, 0.3])
outer_dom  = dde.geometry.CSGDifference(farfield, inner_rec)
outer_dom  = dde.geometry.CSGDifference(outer_dom, cylinder)
inner_dom  = dde.geometry.CSGDifference(inner_rec, cylinder)

farfield_time = dde.geometry.GeometryXTime(farfield, time)
cylinder_time = dde.geometry.GeometryXTime(cylinder, time)
inner_dom_time = dde.geometry.GeometryXTime(inner_dom, time)
outer_dom_time = dde.geometry.GeometryXTime(outer_dom, time)


inner_points = inner_dom_time.random_points(10000)
outer_points = outer_dom_time.random_points(40000)

farfield_points = farfield_time.random_boundary_points(1280)
cylinder_points = cylinder_time.random_boundary_points(250)

# inner_points = inner_dom.random_points(10000)
# outer_points = outer_dom.random_points(40000)

# farfield_points = farfield.random_boundary_points(1280)
# cylinder_points = cylinder.random_boundary_points(250)

points = np.append(inner_points, outer_points, axis = 0)
points = np.append(points, farfield_points, axis = 0)
points = np.append(points, cylinder_points, axis = 0)


# Boundaries definition
def boundary_farfield_inlet(x, on_boundary):
    return on_boundary and np.isclose(x[0], xmin)

def boundary_farfield_top_bottom(x, on_boundary):
    return on_boundary and (np.isclose(x[1], ymax) or np.isclose(x[1], ymin))

def boundary_farfield_outlet(x, on_boundary):
    return on_boundary and np.isclose(x[0], xmax)

def boundary_cylinder(x, on_boundary):
    return on_boundary and cylinder_time.on_boundary(x)

def boundary_init(x, on_inital):
    return on_inital
    
# Boundary values definition
def fun_u_inlet(x, y, _):
    return dde.grad.jacobian(y, x, i = 0, j = 1) - 4 * umax * (ymax - x[:, 1:2]) * x[:, 1:2] / ymax**2

def fun_no_slip_u(x, y, _):
    return dde.grad.jacobian(y, x, i = 0, j = 1)

def fun_no_slip_v(x, y, _):
    return - dde.grad.jacobian(y, x, i = 0, j = 0)

def funP(x):
    return 0.0


# Boundary conditions assembly   
bc_inlet_u = dde.OperatorBC(geom_time, fun_u_inlet, boundary_farfield_inlet)
bc_inlet_v = dde.OperatorBC(geom_time, fun_no_slip_v, boundary_farfield_inlet)

bc_top_bottom_u = dde.OperatorBC(geom_time, fun_no_slip_u, boundary_farfield_top_bottom)
bc_top_bottom_v = dde.OperatorBC(geom_time, fun_no_slip_v, boundary_farfield_top_bottom)

bc_outlet_p = dde.DirichletBC(geom_time, funP, boundary_farfield_outlet, component = 1)

bc_cylinder_u = dde.OperatorBC(geom_time, fun_no_slip_u, boundary_cylinder)
bc_cylinder_v = dde.OperatorBC(geom_time, fun_no_slip_v, boundary_cylinder)

ic_u = dde.OperatorBC(geom_time, fun_u_inlet, boundary_init)
ic_v = dde.OperatorBC(geom_time, fun_no_slip_v, boundary_init)


bcs = [bc_inlet_u, bc_inlet_v, bc_top_bottom_u, bc_top_bottom_v, bc_outlet_p, bc_cylinder_u, bc_cylinder_v, ic_u, ic_v]



# Problem setup
data = dde.data.TimePDE(geom_time, navier_stokes, bcs, num_domain = 0, num_boundary = 0, num_test = 5000, anchors = points)

plt.figure(figsize = (16, 9))
plt.scatter(data.train_x_all[:,0], data.train_x_all[:,1], s = 0.3)
plt.axis('equal')


dde.config.set_default_float('float64')
# Neural network definition
layer_size  = [3] + [40] * 8 + [5]
activation  = 'tanh' 
initializer = 'Glorot uniform'

net = dde.nn.FNN(layer_size, activation, initializer)



# Model definition
model = dde.Model(data, net)
model.compile(optimizer = 'adam', lr = 5e-4, loss_weights = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2]) # Giving more weight to bcs (actually no needed for hard imposed ones)
    
losshistory, train_state = model.train(epochs = 10000, display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)

model.compile(optimizer = 'L-BFGS-B', loss_weights = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2])
model.train_step.optimizer_kwargs = {'options': {'maxcor': 50, 
                                                   'ftol': 1.0 * np.finfo(float).eps, 
                                                   'maxfun':  25000, 
                                                   'maxiter': 25000, 
                                                   'maxls': 50}}
losshistory, train_state = model.train(display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)

time_cylinder_wake.txt

[lululxvi]

There should be something wrong with the code. It is a long code, and you need first figure out which part makes the error.