Skip to content

Commit

Permalink
Adding eigenvector based loss function
Browse files Browse the repository at this point in the history
  • Loading branch information
@jloveric
jloveric committed Feb 12, 2023
1 parent ee341c6 commit eec6f4f
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 38 deletions.
4 changes: 4 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,10 @@ another with contact discontinuity and correct velocity, 8 layers and 10th order
```
python examples/high_order_euler.py mlp.hidden.width=20 max_epochs=10000 mlp.segments=2 mlp.n=2 mlp.hidden.layers=8 factor=0.025 mlp.layer_type=continuous optimizer.patience=200 mlp.input.segments=10 batch_size=2048 form=primitive loss_weight.discontinuity=0.0 loss_weight.interior=1.0e-1 optimizer=adamw mlp.normalize=True mlp.rotations=4 gradient_clip=5.0e-1 loss_weight.boundary=10 loss_weight.initial=10 data_size=10000 mlp.resnet=False refinement.type=p_refine refinement.epochs=500 refinement.target_n=11
```
This one produced a pretty solid rarefaction wave, could be steeper (developed at 6th order) using "waves" approach.
```
python examples/high_order_euler.py mlp.hidden.width=20 max_epochs=10000 mlp.segments=2 mlp.n=2 mlp.hidden.layers=12 factor=0.025 mlp.layer_type=continuous optimizer.patience=200 mlp.input.segments=10 batch_size=2048 form=primitive loss_weight.discontinuity=0.0 loss_weight.interior=1.0e-1 optimizer=adamw mlp.normalize=True mlp.rotations=4 gradient_clip=5.0e-1 loss_weight.boundary=10 loss_weight.initial=10 data_size=10000 mlp.resnet=False refinement.type=p_refine refinement.epochs=500 refinement.target_n=20 refinement.start_n=2
```
## Training
High order MLP
```
Expand Down
140 changes: 106 additions & 34 deletions euler_eigensystem.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"{u - sqrt(p*r*y)/r: 1, u + sqrt(p*r*y)/r: 1, u: 1}\n"
"{-sqrt(c2) + u: 1, sqrt(c2) + u: 1, u: 1}\n"
]
}
],
Expand All @@ -29,13 +29,14 @@
"p=symbols('p')\n",
"y = symbols('y')\n",
"x0=symbols('x0')\n",
"A = Matrix([[u, r, 0], [0, u, 1/r],[0, y*p, u]])\n",
"c2=symbols('c2')\n",
"A = Matrix([[u, r, 0], [0, u, 1/r],[0, y*p, u]]).subs(p,c2*r/y)\n",
"print(A.eigenvals())"
]
},
{
"cell_type": "code",
"execution_count": 26,
"execution_count": 2,
"id": "3d16c5b9",
"metadata": {},
"outputs": [
Expand All @@ -48,21 +49,21 @@
" [1],\n",
" [0],\n",
" [0]])]),\n",
" (u - sqrt(p*r*y)/r,\n",
" (-sqrt(c2) + u,\n",
" 1,\n",
" [Matrix([\n",
" [ r/(p*y)],\n",
" [-u/(p*y) + (u - sqrt(p*r*y)/r)/(p*y)],\n",
" [ 1]])]),\n",
" (u + sqrt(p*r*y)/r,\n",
" [ 1/c2],\n",
" [-u/(c2*r) + (-sqrt(c2) + u)/(c2*r)],\n",
" [ 1]])]),\n",
" (sqrt(c2) + u,\n",
" 1,\n",
" [Matrix([\n",
" [ r/(p*y)],\n",
" [-u/(p*y) + (u + sqrt(p*r*y)/r)/(p*y)],\n",
" [ 1]])])]"
" [ 1/c2],\n",
" [-u/(c2*r) + (sqrt(c2) + u)/(c2*r)],\n",
" [ 1]])])]"
]
},
"execution_count": 26,
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -74,7 +75,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 3,
"id": "e890bbf2",
"metadata": {},
"outputs": [],
Expand All @@ -93,7 +94,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 4,
"id": "babd6973",
"metadata": {},
"outputs": [],
Expand All @@ -103,23 +104,23 @@
},
{
"cell_type": "code",
"execution_count": 32,
"execution_count": 5,
"id": "abba7431",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1\\\\0 & - \\frac{\\sqrt{p r y}}{r^{2}} & \\frac{\\sqrt{p r y}}{r^{2}}\\\\0 & \\frac{p y}{r} & \\frac{p y}{r}\\end{matrix}\\right]$"
"$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1\\\\0 & - \\frac{\\sqrt{c_{2}}}{r} & \\frac{\\sqrt{c_{2}}}{r}\\\\0 & c_{2} & c_{2}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1, 1, 1],\n",
"[0, -sqrt(p*r*y)/r**2, sqrt(p*r*y)/r**2],\n",
"[0, p*y/r, p*y/r]])"
"[1, 1, 1],\n",
"[0, -sqrt(c2)/r, sqrt(c2)/r],\n",
"[0, c2, c2]])"
]
},
"execution_count": 32,
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -130,23 +131,23 @@
},
{
"cell_type": "code",
"execution_count": 31,
"execution_count": 6,
"id": "ebcdc157",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 0 & - \\frac{r}{p y}\\\\0 & - \\frac{r \\sqrt{p r y}}{2 p y} & \\frac{r}{2 p y}\\\\0 & \\frac{r^{2}}{2 \\sqrt{p r y}} & \\frac{r}{2 p y}\\end{matrix}\\right]$"
"$\\displaystyle \\left[\\begin{matrix}1 & 0 & - \\frac{1}{c_{2}}\\\\0 & - \\frac{r}{2 \\sqrt{c_{2}}} & \\frac{1}{2 c_{2}}\\\\0 & \\frac{r}{2 \\sqrt{c_{2}}} & \\frac{1}{2 c_{2}}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1, 0, -r/(p*y)],\n",
"[0, -r*sqrt(p*r*y)/(2*p*y), r/(2*p*y)],\n",
"[0, r**2/(2*sqrt(p*r*y)), r/(2*p*y)]])"
"[1, 0, -1/c2],\n",
"[0, -r/(2*sqrt(c2)), 1/(2*c2)],\n",
"[0, r/(2*sqrt(c2)), 1/(2*c2)]])"
]
},
"execution_count": 31,
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -158,23 +159,23 @@
},
{
"cell_type": "code",
"execution_count": 30,
"execution_count": 7,
"id": "ce40ea67",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}u & 0 & 0\\\\0 & u - \\frac{\\sqrt{p r y}}{r} & 0\\\\0 & 0 & \\frac{p y}{\\sqrt{p r y}} + u\\end{matrix}\\right]$"
"$\\displaystyle \\left[\\begin{matrix}u & 0 & 0\\\\0 & - \\sqrt{c_{2}} + u & 0\\\\0 & 0 & \\sqrt{c_{2}} + u\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[u, 0, 0],\n",
"[0, u - sqrt(p*r*y)/r, 0],\n",
"[0, 0, p*y/sqrt(p*r*y) + u]])"
"[u, 0, 0],\n",
"[0, -sqrt(c2) + u, 0],\n",
"[0, 0, sqrt(c2) + u]])"
]
},
"execution_count": 30,
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -185,10 +186,81 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 8,
"id": "a2f56523",
"metadata": {},
"outputs": [],
"source": [
"v0=symbols('v0')\n",
"v1=symbols('v1')\n",
"v2=symbols('v2')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "b1b73d91",
"metadata": {},
"outputs": [],
"source": [
"v=Matrix([v0,v1,v2])"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "f31fbefa",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}v_{0} - \\frac{v_{2}}{c_{2}}\\\\\\frac{v_{2}}{2 c_{2}} - \\frac{r v_{1}}{2 \\sqrt{c_{2}}}\\\\\\frac{v_{2}}{2 c_{2}} + \\frac{r v_{1}}{2 \\sqrt{c_{2}}}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ v0 - v2/c2],\n",
"[v2/(2*c2) - r*v1/(2*sqrt(c2))],\n",
"[v2/(2*c2) + r*v1/(2*sqrt(c2))]])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"transform = lff*v\n",
"transform"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "bd96fd7c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'Matrix([[v0 - v2/c2], [v2/(2*c2) - r*v1/(2*sqrt(c2))], [v2/(2*c2) + r*v1/(2*sqrt(c2))]])'"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"str(transform)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "431dff6f",
"metadata": {},
"outputs": [],
"source": []
}
],
Expand All @@ -208,7 +280,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
"version": "3.10.7"
}
},
"nbformat": 4,
Expand Down
38 changes: 34 additions & 4 deletions neural_network_pdes/euler.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,22 @@ def __getitem__(self, idx):
return self.input[idx], self.output[idx]


def by_left_eigenvector(rho, p, gamma, v0, v1, v2):
"""
Multiply a vector - in this case the pde, by the
left eigenvector assuming (continuity, velocity, pressure) equations
in that order. This is the primitive left eigenvector as derived
in the jupyer notebook
"""
c2 = torch.clamp(gamma * p / rho, 0.01)

return (
v0 - v2 / c2,
v2 / (2 * c2) - rho * v1 / (2 * torch.sqrt(c2)),
v2 / (2 * c2) + rho * v1 / (2 * torch.sqrt(c2)),
)


def initial_condition_loss(outputs: Tensor, targets: Tensor):
diff = targets - outputs
return torch.dot(diff.flatten(), diff.flatten())
Expand Down Expand Up @@ -92,15 +108,29 @@ def interior_loss(
discontinuity = 1.0 / (eps * (torch.abs(ux) - ux) + 1)

c2 = gamma * p / r
# Note, the equations below are multiplied by r to reduce the loss.

delta = 1e-2

continuity_equation = rt / scale_t + (u * rx + r * ux) / scale_x
velocity_equation = ut / scale_t + (u * ux + (1 / r) * px) / scale_x
pressure_equation = pt / scale_t + (r * c2 * ux + u * px) / scale_x

l1, l2, l3 = by_left_eigenvector(
rho=r,
p=p,
gamma=gamma,
v0=continuity_equation,
v1=velocity_equation,
v2=pressure_equation,
)

r_eq = torch.stack(
[
rt / scale_t + (u * rx + r * ux) / scale_x,
l1,
# Normalize the value before
# (r * ut + r * u * ux + px),
ut / scale_t + (u * ux + (1 / r) * px) / scale_x,
pt / scale_t + (r * c2 * ux + u * px) / scale_x,
l2,
l3,
]
)

Expand Down

0 comments on commit eec6f4f

Please sign in to comment.