lululxvi / deepxde

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

inverse problem #178

Closed Xilun1995 closed 3 years ago

Xilun1995 commented 3 years ago

Hello Lulu, I want to know that is there any possibility we can set variable as a function based on location and time as below image

Xilun1995 commented 3 years ago

I am trying to solve c(r)

lululxvi commented 3 years ago

Yes, you can. See "Q: Define an inverse problem to solve unknown parameters/fields in the PDEs or initial/boundary conditions." at FAQ.

Xilun1995 commented 3 years ago

Thanks for your reply. But I still have a followup question. Since C(r) is a function in terms of r only, but the p is a function in terms of r and t. How can I define c (y[:, 1:2]) as a function only in terms of r in the code.

lululxvi commented 3 years ago

It is doable, but slightly tricky. You can use apply_output_transform to modify the network output, see "Q: By default, initial/boundary conditions are enforced in DeepXDE as soft constraints. How can I enforce them as hard constraints?" for how to use it. Specifically, in apply_output_transform, you manually construct a network whose input is only r, and then use this to replace the original network output as C. See "https://deepxde.readthedocs.io/en/latest/_modules/deepxde/maps/fnn.html#FNN" on how to construct a network.

Xilun1995 commented 3 years ago

Hello lulu, I add the constrains to this problem as below and I attached my code. Since I did not get the correct results, could you please tell me if I am right ? ( I used Cartesian coordinate(x and y) instead of the cylindrical coordinate shown in equation). Thanks in advance. image image image image

geom = dde.geometry.Rectangle([-domain_size, -domain_size], [domain_size, domain_size]) timedomain = dde.geometry.TimeDomain(0, time_size) geomtime = dde.geometry.GeometryXTime(geom, timedomain)

def init_func(x, y): dy_x = tf.gradients(y, x)[0] dy_t = dy_x[:, 2:3] return dy_t

def boundary_initial(x, on_initial): return on_initial

bc = dde.NeumannBC(geomtime,lambda x: 0, lambda _, on_boundary: onboundary) ic = dde.IC( geomtime, lambda x: 0, lambda , on_initial: oninitial ) ic2 = dde.OperatorBC(geomtime, lambda x, y, : init_func(x,y), boundary_initial)

lululxvi commented 3 years ago

The BC and IC look good.

Xilun1995 commented 3 years ago

It is doable, but slightly tricky. You can use apply_output_transform to modify the network output, see "Q: By default, initial/boundary conditions are enforced in DeepXDE as soft constraints. How can I enforce them as hard constraints?" for how to use it. Specifically, in apply_output_transform, you manually construct a network whose input is only r, and then use this to replace the original network output as C. See "https://deepxde.readthedocs.io/en/latest/_modules/deepxde/maps/fnn.html#FNN" on how to construct a network.

I am wondering if I can set C ( sound of speed ) as a matrix( containing 256*256 tf.variable s) instead of using a function. Is there any way to do that ? I really appreciate that.

Xilun1995 commented 3 years ago

It is doable, but slightly tricky. You can use apply_output_transform to modify the network output, see "Q: By default, initial/boundary conditions are enforced in DeepXDE as soft constraints. How can I enforce them as hard constraints?" for how to use it. Specifically, in apply_output_transform, you manually construct a network whose input is only r, and then use this to replace the original network output as C. See "https://deepxde.readthedocs.io/en/latest/_modules/deepxde/maps/fnn.html#FNN" on how to construct a network.

Also, I already know how to add hard constraint to my problem. It is really smart trick. My hard constraint is c(x0,y0,t1) = c(x0,y0,t2). c is always the same, but I did not find a good solution for this constraint, could you please give me some hint?

lululxvi commented 3 years ago
Xilun1995 commented 3 years ago

Thanks. I just follow your suggestion. But the result is really bad not even close. Is there any errors in my apply_output_transform. I attached my code. net.apply_output_transform( lambda x, y : tf.concat( [y[:, 0:1] , tf.layers.dense(tf.layers.dense(tf.layers.dense(tf.layers.dense(tf.layers.dense(x[:,0:2], 10, tf.nn.tanh), 10, tf.nn.tanh), 30, tf.nn.tanh), 30, tf.nn.tanh), 1, tf.nn.tanh) ], axis = 1 ) )

lululxvi commented 3 years ago

Should the last layer use activation?

Xilun1995 commented 3 years ago

THANKS!! One more question: I run the tests several times but Every time I got the different results and one of them is correct. Is there any way to improve that. Also, for the 2d wave equation with time, do you think my number of training points and size of NN is correct or not?

num_domain= 12000, num_boundary=1000, num_initial= 2000

[3] + [40] * 7 + [1]

lululxvi commented 3 years ago
Xilun1995 commented 3 years ago

Hello Lulu, when I used 6000 training points, the networks works. However, when I increase the number of training points to 60000, I get the following errors. Can you help me to debug? image

lululxvi commented 3 years ago

It is out of memory.

AlirezaPa commented 2 years ago

Thank you @lululxvi for your library. It is very helpful for our research. I am writing a simple inverse code for 1D-wave equation. First, I modified my code, which I provide below, based on the discussion #114 . Currently, C in my code is a function of x and t, but as we know C is just function of x, and my issue is the same as @Xilun1995 here. I read all above discussions and API, but I still do not know how to write "apply_output_transform" and "tf.layers.dense". Would you please take a look to my code and give me some advises that how to modify my code based on C(x)? I appreciate for your time.

Inverse Wave equation with Dirichlet Boundary Conditions

import deepxde as dde import numpy as np

A = 1 Lright= 1 Lleft = -1 Ldomain = Lright - Lleft Ttotal = 1

def pde(x, y): u, C = y[:, 0:1], y[:, 1:2] du_tt = dde.grad.hessian(u, x, i=1, j=1) du_x = dde.grad.jacobian(u, x, j=0) du_xx = dde.grad.jacobian((C * 2)du_x, x, j=0)
return du_tt - du_xx

def func(x): return np.sin(np.pi A x[:, 0:1]) np.cos(np.pi A 1 x[:, 1:])

def boundary_r(x, on_boundary): return on_boundary and np.isclose(x[0], Lright)

def boundary_l(x, on_boundary): return on_boundary and np.isclose(x[0], Lleft)

geom = dde.geometry.Interval(Lleft, Lright) timedomain = dde.geometry.TimeDomain(0, Ttotal) geomtime = dde.geometry.GeometryXTime(geom, timedomain)

bc_1 = dde.DirichletBC(geomtime, lambda x: 0, boundary_r, component=0) bc_2 = dde.DirichletBC(geomtime, lambda x: 0, boundary_l, component=0) ic1 = dde.IC(geomtime, lambda x: np.sin(np.pi A x[:, 0:1]), lambda , on_initial: on_initial, component=0) ic2 = dde.OperatorBC(geomtime, lambda x, u, : dde.grad.jacobian(u, x, i=0, j=1), lambda x, _: np.isclose(x[1], 0))

observe_x = np.vstack((np.full((300), 0.75), np.linspace(-1, 1, num=300))).T observe_u = dde.icbc.PointSetBC(observe_x, func(observe_x), component=0)

data = dde.data.TimePDE( geomtime, pde, [bc_1, bc_2, ic_1, ic_2, observe_u], num_domain=1000, num_boundary=300, num_initial=300, anchors=observe_x, num_test=10000, )

net = dde.nn.PFNN([2] + [50, 50] * 3 + [2], "tanh", "Glorot normal") model = dde.Model(data, net)

model.compile("adam", lr=0.0001, loss_weights=[0.1, 10, 10, 10, 10, 10]) losshistory, train_state = model.train(epochs=10000) dde.saveplot(losshistory, train_state, issave=True, isplot=True)

lululxvi commented 2 years ago
* I mean you construct a new network in `apply_output_transform` as
def apply_output_transform(inputs, outputs):
    p = outputs[:, 0:1]
    x = inputs[:, 0:1]
    C = FNN(x)
    return tf.concat([p, C], axis=1)

In C = FNN(x), construct a network using tf.layers.dense.

This is the answer. You can also find examples at FAQ "Q: A standard network doesn’t work for me, and I want to implement a network with a special structure/property."

hannanmustajab commented 11 months ago

@AlirezaPa Hey, I am also solving the same problem. How did you manage to implement this in the end ? Can you please share your code snippet for that ? Also, did you use tf backend only ?

I also have spatially varying velocity model.