Whole domain heat flux estimation as function of time

[rodrigogdourado]

Hello,

I'm trying to estimate the heat flux boundary condition given by the discrete values q = [q_1,q_2, ..., q_N], in the whole time domain t = [t_1, t_2, ... , t_N] using known data Y = [Y_1, Y_2, ..., Y_N] at a point within the domain. I've already tried several ways without success. Now, I'm defining q initial guess as a vector with N components: q = dde.Variable(np.ones[(len(Y),1], dtype='float32')), but I'm not sure if it works since it doesn't seem like the algorithm update the q components over the iterations.

I would like to know if someone has an idea how to solve this problem.

Here is the problem (estimate q using Y):

Here is the code:

import deepxde as dde
# import tensorflow as tf
from deepxde.backend import tf
import numpy as np
import numbers

#%% parameters
# Y = np.dtype('float32')
# q = np.dtype('float32')
dados = np.loadtxt('C:\DOURADO\PYTHON_TREINAMENTO\Y_whole.txt')
Y = dados.reshape(-1,1)

q = dde.Variable(np.ones([len(Y),1],dtype='float32'))

L = 1;
time =1;
T0 = 0;
k = 1.2;
rho = 1;
cp = 1;
alpha = k/(rho*cp)
dt = 0.05

#anchors
# tt = np.arange(0,time+0.001,time/50).reshape(-1,1)
tt =np.arange(0,time+0.01,dt).reshape(-1,1)
xx = np.zeros([len(tt),1])
ancora_x = np.hstack([xx,tt])   #anchor for heat flux boundary condition
x_sensor = np.ones([len(Y),1]).reshape(-1,1)*L
observe_x = np.hstack([x_sensor,tt])    #anchor for sensor location (inverse problem)

#heat flux boundary condition using PointSetBCDot

class PointSetBCDot(object):

    def __init__(self, points, values, component, j):
        self.points = np.array(points)
        if not isinstance(values, numbers.Number) and values.shape[1] != 1:
            raise RuntimeError(
                "PointSetBC should output 1D values. Use argument 'component' for different components."
            )
        self.j = j
        self.values = values
        self.component = component

    def first_normal_derivative(self, X, inputs, outputs, beg, end):
        dydt = dde.grad.jacobian(outputs, inputs, i=self.component, j=self.j)[beg:end]
        return tf.reduce_sum(dydt, axis=1, keepdims=True)

    def collocation_points(self, X):
        return self.points

    def error(self, X, inputs, outputs, beg, end):
        #return self.first_normal_derivative(X, inputs, outputs, beg, end) - tf.constant(self.values)
        return self.first_normal_derivative(X, inputs, outputs, beg, end) + self.values/k

#domain

geom = dde.geometry.Interval(0, L)
timedomain = dde.geometry.TimeDomain(0, time)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

#PDE and loss

def pde(x,y):
    dy_xx = dde.grad.hessian(y, x, i=0,j=0)
    dy_t = dde.grad.jacobian(y, x, i=0,j=1)
    return dy_t-alpha*dy_xx

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

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

#Boundary

bc_l = PointSetBCDot(ancora_x,q,component=0,j=0)    #heat flux boundary condition
bc_r = dde.icbc.NeumannBC(geomtime, lambda x:0, boundary_r)
ic = dde.icbc.IC(geomtime, lambda x:T0, lambda _, on_initial:on_initial)
observe_y = dde.icbc.PointSetBC(observe_x, Y, component=0)  # reference for inverse problem

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc_l,bc_r,ic,observe_y],
    num_domain=500,
    num_boundary=200,
    num_initial=80,
    anchors = observe_x,
    num_test=None,
    )

layer_size = [2]+[64]*3+[1]
net = dde.nn.FNN(layer_size,"tanh","Glorot normal")
#net = dde.nn.FNN(layer_size,"tanh","Glorot normal")
model = dde.Model(data,net)

model.compile("adam",lr=0.1e-3,loss_weights=[1,2,1,1,10],external_trainable_variables=q)
variable = dde.callbacks.VariableValue(q, period=1000)
losshistory, train_state = model.train(iterations=40000)

[lululxvi]

#987