# Variational loss of hp-Variational Physics Informed Neural Networks for 2D-Poisson Equation in Tensorflow

I am trying to reproduce the results from the hp-VPINN paper (https://arxiv.org/pdf/2003.05385.pdf) on tensorflow (v1) for Poisson's equation, particularly the two-dimensional Poisson equation.

In one dimension the variational loss is given by $$L^p = \sum_{e=1}^{N_{el}}\frac{1}{K^{(e)}}\sum_{k=1}^{K^{(e)}}\left|\mathcal{R}_k^{(e)}\right|^2$$

where $$N_{el}$$ is the number of elements and $$K^{(e)}$$ is the number of test functions used (in this case corresponding to different degrees of Legendre polynomials)

where $$\mathcal{R}_k^{(e)}$$ is $$\mathcal{R}_k^{(e)} = \int_{x_{e-1}}^{x_e} \frac{du_{NN}(x)}{dx} \frac{dv_k^{(e)}(x)}{dx}dx - \left. \frac{du_{NN}(x)}{dx}v_k^{(e)}(x)\right|_{x_{e-1}}^{x_{e}} - F_k^{(e)}$$

with $$u_{NN}$$ being the network's output and and $$v_k^{(e)}$$ the test function. The boundary term vanishes since the test functions are such that they are zero on the element boundaries and $$F_k^{(e)} = \int_{x_{e-1}}^{x_e} f(x)v_k^{(e)}dx$$.

To calculate it Gauss-Lobatto quadrature is used giving the following code

import tensorflow as tf
import numpy as np
u_nn_element = tf.reshape(
tf.stack(
[
tf.reduce_sum(
)
for i in range(n_test_functions)
]
),
(-1, 1),
)
f_element = jacobian * np.asarray(
)
residual_element = u_nn_element - f_element
loss_element = tf.reduce_mean(tf.square(residual_nn_element))
total_varloss += loss_element


When combined with the boundary loss (mean squared error for the boundary points at $$[-1, 1]$$)

In two dimensions the variational loss is given by

$$L= \sum_{e_x=1}^{N_{el_x}}\sum_{e_y=1}^{N_{el_y}}\frac{1}{K_1K_2}\sum_{k=1}^{K_1K_2}\left|\mathcal{R}_k^{(e_xe_y)}\right|^2$$

with a residual given by

$${}^{(2)} \mathcal{R}_{k_1k_2}^{(e_xe_y)} = - \int_{x_{e-1}}^{x_e}\int_{y_{e-1}}^{y_e} \left(\frac{\partial u_{NN}(x)}{\partial x} \frac{d\phi_{k_1}^{(e_x)}(x)}{dx}\phi_{k_2}^{(e_y)}(y) + \frac{\partial u_{NN}}{\partial y}\phi_{k_1}^{(e_x)}(x)\frac{d\phi_{k_2}^{(e_y)}(y)}{dy} \right) dxdy \\ - F_k^{(e)}$$

with the test functions $$\phi_{k_1} \phi_{k2}$$

My code attempt has been

l = [
tf.reduce_sum(
* (
+ du_nn_dy
)
)
for i in range(n_test_functions_x)
for j in range(n_test_functions_y)
]
u_nn_element = tf.reshape(
tf.stack(l),
(-1, 1),
)

# element boundaries
jacobian_x = (x1 - x0) / 2
jacobian_y = (y1 - y0) / 2

f = [
for t1 in self.test_functions
for t2 in self.test_functions
]
self.f = jacobian_x * jacobian_y * np.asarray(f)
self.f = self.f.flatten()
self.f = self.f[:, None]

residual_nn_element = u_nn_element - f_element
loss_element = tf.reduce_mean(tf.square(residual_nn_element))
self.varloss_total += loss_element


though my variational error doesn't converge and the net only learns the boundary values. How is the double integral supposed to look in code?