lululxvi / deepxde

A library for scientific machine learning and physics-informed learning
https://deepxde.readthedocs.io
GNU Lesser General Public License v2.1
2.64k stars 739 forks source link

BC Over Geometry Part #59

Closed sumantkrsoni closed 2 years ago

sumantkrsoni commented 4 years ago

hii @lululxvi , I am solving Navier stoke equation which three systems of PDE equation:

while compiling it shows an error for concatenation:

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 2000 and the array at index 1 has size 1000

Also, in the Geometry part, I have taken a circle with radius 0.2 centered at the center of rectangle:

Whose code is:

    geom1 = dde.geometry.Rectangle([0,0], [ 12, 5 ])
    geom2 = dde.geometry.Disk( [4.0 , 2.5], 0.2 )

    geom = geom1-geom2

I want to put the values of variables over the circle boundary part : u = 0; (component = 0) v = 0; ( component =1)

I have coded it as:

bc_circle_u = dde.DirichletBC(geom2, lambda x: np.zeros( (len(x) ,1 ) ), lambda _, on_boundary: on_boundary, component =0)

bc_circle_v = dde.DirichletBC(geom2, lambda x: np.zeros( (len(x) ,1 ) ), lambda _, on_boundary: on_boundary, component =1)
lululxvi commented 4 years ago

Check #49 to see how to define the BC

sumantkrsoni commented 4 years ago

Thanks @lululxvi ,

It helped me to figure out my doubts. But, still a couple of doubt I have

I have compiled my code for Navier Stokes Equation:

also it shows an error as:

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 2000 and the array at index 1 has size 1000

Kindly suggest


from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np

import deepxde as dde

import tensorflow as tf

def main():

    def pde( x, y ):

        eps = 0.00000001
        Re = 600.0 
        nu = 1/Re

        u, v, p = y[:, 0:1], y[:, 1:2 ], y[:, 2:3]

        du_x = tf.gradients( u , x )[0]
        dv_x = tf.gradients( v , x )[0]
        dp_x = tf.gradients( p , x )[0]

        dp_x, dp_y = dp_x[:, 0:1 ], dp_x[:, 1:2 ]

        du_x, du_y , du_t = du_x[:, 0:1], du_x[:, 1:2 ], du_x[:, 2:3 ]

        dv_x, dv_y, dv_t = dv_x[:, 0:1], dv_x[ :, 1:2 ], dv_x[:, 2:3 ]

        du_xx = tf.gradients( du_x, x )[0][:, 0:1 ]
        du_yy = tf.gradients( du_y, x )[0][:, 1:2 ]

        dv_xx = tf.gradients( dv_x, x)[0][:, 0:1 ]
        dv_yy = tf.gradients( dv_y, x)[0][:, 1:2 ]

        return [ du_x + du_y + eps * p, 
                du_t + u * du_x + v * du_y - nu * (du_xx + du_yy ) + dp_x ,
                dv_t + u * dv_x + v * dv_y - nu * ( dv_xx + dv_yy ) + dp_y ]

    # function definition

    def func(x):

        return np.zeros((len(x),1))

    def func_g(x):

        Ly = 5.0
        return 4 * ( x[ : , 1:2 ] / Ly )* ( 1 - x[ : , 1:2]/ Ly ) 

    #return np.ones((len(x), 1 ))

    rect_geom = dde.geometry.Rectangle([0,0], [ 12, 5 ])
    circle_geom = dde.geometry.Disk( [4.0 , 2.5], 0.2 )

    domain_geom = dde.geometry.CSGDifference( rect_geom, circle_geom )

    timedomain = dde.geometry.TimeDomain( 0, 5 )

    geomtime = dde.geometry.GeometryXTime( domain_geom , timedomain)

    # boundary definition of a unit square domain

    def boundary_left( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x) and np.isclose(x[0], 0 )

    def boundary_right( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x) and np.isclose(x[0], 12 )

    def boundary_upper( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x) and np.isclose(x[1], 5 )

    def boundary_lower( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x) and np.isclose(x[1], 0 )

    def boundary_circle( x, on_boundary ):

        return on_boundary and circle_geom.on_boundary( x ) 

    # defining BC for u, v, and p (component = 0, 1, 2 resp assign)

    bc_left_u = dde.DirichletBC( geomtime, func_g , boundary_left, component = 0 )
    bc_left_v = dde.DirichletBC( geomtime, func , boundary_left, component = 1 )

    bc_right_u = dde.NeumannBC( geomtime, func, boundary_right, component = 0 )
    bc_right_v = dde.NeumannBC( geomtime, func, boundary_right, component = 1 )
    bc_right_p = dde.DirichletBC(geomtime, func, boundary_right, component = 2 )

    bc_lower_u = dde.DirichletBC( geomtime, func, boundary_lower, component = 0 )
    bc_lower_v = dde.DirichletBC( geomtime, func, boundary_lower, component = 1 )

    bc_upper_u = dde.DirichletBC(geomtime, func, boundary_upper, component = 0 )
    bc_upper_v = dde.DirichletBC(geomtime, func, boundary_upper, component = 1 )

    bc_circle_u = dde.DirichletBC(geomtime, func, boundary_circle, component =0 )
    bc_circle_v = dde.DirichletBC(geomtime, func, boundary_circle, component =1 )

    ic_u = dde.IC(  geomtime, lambda x: np.zeros((len(x),1)) , lambda _, on_initial: on_initial, component = 0 )
    ic_v = dde.IC(  geomtime, lambda x: np.zeros((len(x),1)) , lambda _, on_initial: on_initial, component = 1 )

    #bc_lwr_p = dde.DirichletBC(geom, lambda x: 0.00001 * np.ones((len(x), 1)), boundary_upr, component = 2)

    bc_ic = [bc_left_u, bc_left_v, 
             bc_right_u, bc_right_v, bc_right_p, 
             bc_lower_u, bc_lower_v, 
             bc_upper_u, bc_upper_v, 
             bc_circle_u, bc_circle_v, 
             ic_u, ic_v ]

    data = dde.data.TimePDE(
                     geomtime,
                     pde,
                     bc_ic,
                     num_domain = 10**4,
                     num_boundary = 10**3,
                     num_initial = 8**3,
                     num_test = 20**2
                   )

    layer_size = [3] + [50] * 2 + [3]
    activation = "tanh"
    initializer = "Glorot uniform"

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

    model = dde.Model(data, net)

    model.compile( "adam", lr = 0.001 )

    losshistory, train_state = model.train(epochs = 15000)

    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

if __name__ == "__main__":
    main()
lululxvi commented 4 years ago

The code looks good. Could you first check which IC or BC induces this error?

sumantkrsoni commented 4 years ago

I think BC induces this error.

Trying to find the exact BC which causes this error.

Can you suggest me why this dimensional agreement does not match with the variables?

I have another complex non-dimensional equation for Navier Stokes Equation and for my better understanding I have tried this problem over DeepXDE.

Hopefully, You will help me.

Kindly help @lululxvi @tangqi @smao-astro

sumantkrsoni commented 4 years ago
                 After replacing the collocation point over domain and Boundary: 

I got error as: ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1600 and the array at index 1 has size 800

Different training points:

data = dde.data.TimePDE(
                     geomtime,
                     pde,
                     bc_ic,
                     num_domain = 100,
                     num_boundary = 800,
                     num_initial = 300,
                     num_test = 2000
                   ) 

error is related to the domain_boundary = 800

lululxvi commented 4 years ago

There is a bug. I have fixed the bug. You can download the updated code.

When defining the boundary, using only x and y, but not t, i.e., x[:2], see the example of boundary_circle:

def boundary_circle(x, on_boundary):
    return on_boundary and circle_geom.on_boundary(x[:2])
sumantkrsoni commented 4 years ago

Thanks @lululxvi

Thanks for your support.

I am not getting the code: " You can download the updated code."

Since, x, y & t are using over the domain. (rectangle - circle)

then why should we need to exclude the '' t " variable.

After, updating the code according to your suggestion the same error I am getting:

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1600 and the array at index 1 has size 800

Kindly, suggest .

Thanks.

lululxvi commented 4 years ago
sumantkrsoni commented 4 years ago

Thanks @lululxvi

Further, I got an error in the line

sys.path.insert( 0 , ......)

runfile('D:/PhD_Box/Code/Project_1/Stoke_2_DeepXDE/navierTransient2.py', wdir='D:/PhD_Box/Code/Project_1/Stoke_2_DeepXDE')
  File "D:\PhD_Box\Code\Project_1\Stoke_2_DeepXDE\navierTransient2.py", line 10
    sys.path.insert( 0, "C:\Users\SUMANT\Anaconda3\envs\tensorflow1\Lib\site-packages\deepxde")
                       ^
SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 2-3: truncated \UXXXXXXXX escape

Then, I modified the path with including 'r' before path as:

sys.path.insert( 0, r"C:\Users\SUMANT\Anaconda3\envs\tensorflow1\Lib\site-packages\deepxde")

also with:

sys.path.insert( 0, "C:\\Users\\SUMANT\\Anaconda3\\envs\tensorflow1\\Lib\site-packages\\deepxde") then, that Unicode error get fixed.

But, unfortunately, after compilation, I am getting the same error:

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1600 and the array at index 1 has size 800

Final updated code is given as ( After modifying the different boundary):


from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import sys
sys.path.insert( 0,  r"C:\Users\SUMANT\Anaconda3\envs\tensorflow1\Lib\site-packages\deepxde")
#sys.path.insert( 0,  "C:\\Users\\SUMANT\\Anaconda3\\envs\\tensorflow1\\Lib\\site-packages\\deepxde")

import deepxde as dde
import numpy as np

import tensorflow as tf

def main():

    def pde( x, y ):

        eps = 0.00000001
        Re = 600.0 
        nu = 1/Re

        u, v, p = y[:, 0:1], y[:, 1:2 ], y[:, 2:3]

        du_x = tf.gradients( u , x )[0]
        dv_x = tf.gradients( v , x )[0]
        dp_x = tf.gradients( p , x )[0]

        dp_x, dp_y = dp_x[:, 0:1 ], dp_x[:, 1:2 ]

        du_x, du_y , du_t = du_x[:, 0:1], du_x[:, 1:2 ], du_x[:, 2:3 ]

        dv_x, dv_y, dv_t = dv_x[:, 0:1], dv_x[ :, 1:2 ], dv_x[:, 2:3 ]

        du_xx = tf.gradients( du_x, x )[0][:, 0:1 ]
        du_yy = tf.gradients( du_y, x )[0][:, 1:2 ]

        dv_xx = tf.gradients( dv_x, x)[0][:, 0:1 ]
        dv_yy = tf.gradients( dv_y, x)[0][:, 1:2 ]

        return [ du_x + du_y + eps * p, 
                du_t + u * du_x + v * du_y - nu * (du_xx + du_yy ) + dp_x ,
                dv_t + u * dv_x + v * dv_y - nu * ( dv_xx + dv_yy ) + dp_y ]

    # function definition

    def func(x):

        return np.zeros((len(x),1))

    def func_g(x):

        Ly = 5.0
        return 4 * ( x[ : , 1:2 ] / Ly )* ( 1 - x[ : , 1:2]/ Ly ) 

    #return np.ones((len(x), 1 ))

    rect_geom = dde.geometry.Rectangle([0,0], [ 12, 5 ])
    circle_geom = dde.geometry.Disk( [4.0 , 2.5], 0.2 )

    domain_geom = dde.geometry.CSGDifference( rect_geom, circle_geom )

    timedomain = dde.geometry.TimeDomain( 0, 5 )

    geomtime = dde.geometry.GeometryXTime( domain_geom , timedomain)

    # boundary definition of a unit square domain

    def boundary_left( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x[:2]) and np.isclose(x[0], 0 )

    def boundary_right( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x[:2]) and np.isclose(x[0], 12 )

    def boundary_upper( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x[:2]) and np.isclose(x[1], 5 )

    def boundary_lower( x, on_boundary ):

        return on_boundary and rect_geom.on_boundary(x[:2]) and np.isclose(x[1], 0 )

    def boundary_circle( x, on_boundary ):

        return on_boundary and circle_geom.on_boundary( x[:2] ) 

    # defining BC for u, v, and p (component = 0, 1, 2 resp assign)

    bc_left_u = dde.DirichletBC( geomtime, func_g , boundary_left, component = 0 )
    bc_left_v = dde.DirichletBC( geomtime, func , boundary_left, component = 1 )

    bc_right_u = dde.NeumannBC( geomtime, func, boundary_right, component = 0 )
    bc_right_v = dde.NeumannBC( geomtime, func, boundary_right, component = 1 )
    bc_right_p = dde.DirichletBC( geomtime, func, boundary_right, component = 2 )

    bc_lower_u = dde.DirichletBC( geomtime, func, boundary_lower, component = 0 )
    bc_lower_v = dde.DirichletBC( geomtime, func, boundary_lower, component = 1 )

    bc_upper_u = dde.DirichletBC( geomtime, func, boundary_upper, component = 0 )
    bc_upper_v = dde.DirichletBC( geomtime, func, boundary_upper, component = 1 )

    bc_circle_u = dde.DirichletBC( geomtime, func, boundary_circle, component =0 )
    bc_circle_v = dde.DirichletBC( geomtime, func, boundary_circle, component =1 )

    ic_u = dde.IC(  geomtime, lambda x: np.zeros((len(x),1)) , lambda _, on_initial: on_initial, component = 0 )
    ic_v = dde.IC(  geomtime, lambda x: np.zeros((len(x),1)) , lambda _, on_initial: on_initial, component = 1 )

    #bc_lwr_p = dde.DirichletBC(geom, lambda x: 0.00001 * np.ones((len(x), 1)), boundary_upr, component = 2)

    bc_ic = [ bc_left_u, bc_left_v, 
             bc_right_u, bc_right_v, bc_right_p, 
             bc_lower_u, bc_lower_v, 
             bc_upper_u, bc_upper_v, 
             bc_circle_u, bc_circle_v, 
             ic_u, ic_v ]

    data = dde.data.TimePDE(
                     geomtime,
                     pde,
                     bc_ic,
                     num_domain = 100,
                     num_boundary = 800,
                     num_initial = 300,
                     num_test = 2000
                   )

    layer_size = [3] + [50] * 2 + [3]
    activation = "tanh"
    initializer = "Glorot uniform"

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

    model = dde.Model(data, net)

    model.compile( "adam", lr = 0.001 )

    losshistory, train_state = model.train(epochs = 15000)

    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

if __name__ == "__main__":
    main()
lululxvi commented 4 years ago

Install the new version v0.8.1

sumantkrsoni commented 4 years ago

Thanks @lululxvi -Fortunately, now It works well.

I am trying to predict the values and save its values over domain for different time = 0 and 3 .

and I have coded it as:

  x_geom = domain_geom.uniform_points( 10000, True )
    t_geom = np.ones_like( ( len(x_geom), 1 ) )

    X = np.vstack( ( x_geom, t_geom ) )
    y_geom_output = model.predict( X , operator = pde )

    np.savetxt( "test_predict_data.dat", np.hstack( (X, y_geom_output) ) )

But, I dimensional conflicts occurs here.

ValueError: Cannot feed value of shape (10001, 2) for Tensor 'Placeholder_390:0', which has shape '(?, 3)'

kindly, suggest me to predict the values over domain.

lululxvi commented 4 years ago
sumantkrsoni commented 4 years ago

I have updated the predicted part of the code:

    # uniform points over geometry

    X_geom_uniform = geomtime.uniform_points( 30000, boundary = True )

    #t_geom = np.zeros_like( ( len(x_geom), 1 ))

    #X = np.vstack( ( x_geom, t_geom ) )

    print("Shape of X uniform over geometry is: ", X_geom_uniform.shape)

    y_geom_output = model.predict( X_geom_uniform )

    np.savetxt( "test_predict_uniform.dat", np.hstack( (X_geom_uniform, y_geom_output) ) )

    # random points over geometry
    X_geom_random = geomtime.random_points( 20000, random = "pseudo")

    print("Shape of X rendom over geometry is: ", X_geom_random.shape)

    y_geom_random = model.predict( X_geom_random )

    np.savetxt( "test_predict_random.dat", np.hstack( (X_geom_random, y_geom_random) ) )

Now, almost everything is working fine, except

say, in my case time domain is = (0, 20 ),

then I am trying to plot the contour at time t = 0, 1 , 2, ..., 20.

I have got the test data arranged into x, y, t, u, v, p variable (columnwise).

can, I extract the data at every time step from test.dat file.

the points over domain not as expected, it has less number of points over domain and boundary: cir_rect_geom

the points over domain and boundary are fine: random_points

But, the predicted output value from test_predict_random.dat is totally unordered output for the time.

Kindly, suggest me. @lululxvi

Test_file.zip

lululxvi commented 4 years ago

In this geometry (rectangle - circle), the points are sampled randomly, because it is not clear what the uniform or grid means here. If you want to control the locations of points, you can construct the points manually, instead of using geomtime.uniform_points() or geomtime.random_points(). But rect_geom.uniform_points() could be helpful.

sumantkrsoni commented 4 years ago

I am trying to modify this code for

I am in the process to solve each part one by one.

Thanks

lululxvi commented 4 years ago
sumantkrsoni commented 4 years ago

Thanks @lululxvi

If any model is already trained for epochs = 10000 (say) and further trying to trained further till 15000 epochs then.

According to #57 both codes should put together simultaneously ( training and restoring )

https://github.com/lululxvi/deepxde/issues/57#issuecomment-637228551

This is crucial information that would save my time for the large coupled problem.

Kindly suggest.

smurf-17 commented 3 years ago

@SanjuSoni I tried the flow over cylinder for 100 Reynolds number but while plotting the result I am getting the below plot. Screenshot (616) Could you help me with getting the plot. I am trying to solve this question with Deepxde :- https://drive.google.com/file/d/1FryjovP2g_a_dvl_KHPkNQ9T0jEBpddf/view and I want to capture the von karman shedding

sumantkrsoni commented 3 years ago

@smurf-17 May I know, which code did you compile?

smurf-17 commented 3 years ago

@SanjuSoni Sir this is the link for the code.

smurf-17 commented 3 years ago

@SanjuSoni Sir I am very new to Deepxde so there can be issue with the code too, So it will really be great if you can give me suggestion about the code too.

sumantkrsoni commented 3 years ago

Hii @smurf-17 Your code seems fine for now except few scripts that are mentioned below:

x-component velocity is mismatching from the attached ppt ( please cross-check it).

def func_u(x):
  '''function of u'''
  return x[:,0:1] *0+0.0146

no-need to add an extra term in continuity eqn. (eps*p)

continuity = u_x + u_y + eps*p

is it mandatory to provide the boundary condition inside the circle?

  bc_inside_x= deepxde.DirichletBC(geomtime,func_v,circle_inside,component=0)
  bc_inside_y= deepxde.DirichletBC(geomtime,func_v,circle_inside,component=1)

If the above-mentioned points need to be corrected then, please modify them and re-run your code.

Also, I would like to suggest you

smurf-17 commented 3 years ago

@SanjuSoni Thank you sir I will try it again.

smurf-17 commented 3 years ago

@SanjuSoni sir I corrected the code and tried for higher epochs. But the result is:- Screenshot (630) Sir, if we are doing CSG Difference there should be no training points inside the disk but according to the plot there are training points inside the disk. I even tried plotting in matlab , got almost the same results. Also I got a good loss so I think there is no issue of convergence Best model at step 17000: train loss: 1.11e-07 test loss: 3.28e-08 test metric: []

sumantkrsoni commented 3 years ago

can you please show the training points over the entire computational domain?

smurf-17 commented 3 years ago

@SanjuSoni Sir I am very new to using this, I have no idea how to do that. If possible could you please share the code snippet and I will try to run it.

lululxvi commented 3 years ago

@smurf-17 Your rectangle is very large [0.0,0.0]x[31.5, 25.0], but the inner disk is tiny with radius 0.5. That is why you cannot see the disk in the plot. Try a smaller rectangle.

sumantkrsoni commented 3 years ago

@smurf-17 increase the radius from 0.5 to a significantly large value and monitor the plot. Also, upload the test file in text or dat format.

smurf-17 commented 3 years ago

@SanjuSoni @lululxvi I tried to make the size of the disk and rectangle and got the results but while plotting results on matlab using contourf I am still not getting the disk but when I tried scatter plot the disc was visible. test.xlsx

This is the result while plotting using matplotlib. Could you please tell me how can I get the same plot on Matlab Screenshot (685)

sumantkrsoni commented 3 years ago

Hello @smurf-17 Take a look at the attached rar file and try to modify it according to your problem.

Plotting with the scatter plot in matplotlib, I am able to get a circle inside the rectangle shape.

Also, it would definitely display in Matlab too.

Also, if you need further support feel free to comment here. cyl11.zip

cyl1

smurf-17 commented 3 years ago

@SanjuSoni It gives a very good result on MATLAB as well when I use scatter plot but is it possible to get the contour plot and also is it possible to get plots at different time steps

sumantkrsoni commented 3 years ago

for a prediction of the result at a particular time, you can use "model.predict" to get your answer.

Also, for your contour plot, you may use the data scattered at the domain and the result at those points in python. I have not used Matlab so I can't help you with that.

Hopefully, I could answer your query, however, Dr @lululxvi Sir can answer you more appropriately.

lululxvi commented 3 years ago

@SanjuSoni Thanks for your nice answers. @smurf-17 Now, there is nothing wrong in DeepXDE side. The issue is that the contour plot function can display what you need. I didn't plot this type of thing before, and cannot provide you a quick answer.

smurf-17 commented 3 years ago

@lululxvi okk sir, Thank you, I will try to get it

Karankaren commented 2 years ago

Hi LuLu please can you explain to me,

  1. How can I write the code of continuity condition that you used in your presentation, Outer circle: φ1 = φ2, \sigma_1 ∂φ1/∂n = \sigma_2 ∂φ2/∂n

Inner circle: φ2 = φ3, \sigma_2 ∂φ2/∂n = \sigma_3 ∂φ3/∂n

  1. how can I write the Dirichlet and Nuemann condition in two dimension rectangular domain (Rectangle(xmin=[0, 0], xmax=[1, 1])) with respect reference function i.e. for example if reference function is u(x,y)= cos(x)*cosh(y) then how I write DeepXDE code for Nuemann condition ∂u/∂n on a patr of boundary and Dirichlet condition on other part.
lululxvi commented 2 years ago

Check FAQ.

zyz-rise commented 2 years ago

Thanks @lululxvi

  • Early stopping is working well.
  • Only thing is the saving and restoring part is making me little trouble for me. I am confused at one point.

If any model is already trained for epochs = 10000 (say) and further trying to trained further till 15000 epochs then.

According to #57 both codes should put together simultaneously ( training and restoring )

#57 (comment)

This is crucial information that would save my time for the large coupled problem.

Kindly suggest.

Sir, @SanjuSoni I wonder if you can observe an obvious vortex street in your problem of flow around a cylinder? And could you provide the final code If it is convenient ? thanks a lot

sumantkrsoni commented 2 years ago

Hello @zyz-rise, I am not in touch with this topic for a long time but I could give you a complete code script here.

Apart from the code script, if you need any help then please let me know.

Thanks.

navier_transient.zip