Meshes

This guide demonstrates the use of meshes in continuiti.

In this example, we will load a Gmsh file and train a physics-informed neural operator.

import torch
import pathlib
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation
from continuiti.data import OperatorDataset
from continuiti.data.mesh import Gmsh
from continuiti.operators import DeepONet
from continuiti.trainer import Trainer
from continuiti.pde import div, grad, PhysicsInformedLoss

We use the Gmsh class to read a mesh from a file generated by Gmsh.

meshes_dir = pathlib.Path.cwd().joinpath("..", "data", "meshes")
gmsh_file = meshes_dir.joinpath("mediterranean.msh").as_posix()

mesh = Gmsh(gmsh_file)
vertices = mesh.get_vertices()

# We use longitude and latitude only and scale them to be in the range [-1, 1]
vertices = vertices[1:, :]
vertices -= vertices.mean(dim=1, keepdim=True)
vertices /= vertices[0, :].abs().max()

Problem Statement

Let's assume we have a set of temperature measurements and want to reconstruct a physically meaningful temperature distribution in the whole mediterranean sea. We define a data set by generating random sensor inputs and plot the first entry.

num_observation = 1
num_sensors = 5

x = torch.zeros(num_observation, 2, num_sensors)
u = torch.zeros(num_observation, 1, num_sensors)
y = torch.zeros(num_observation, 2, num_sensors)
v = torch.zeros(num_observation, 1, num_sensors)

for i in range(num_observation):

    # Select random sensor positions
    idx = torch.randperm(vertices.shape[1])[:num_sensors]
    x[i] = vertices[:, idx]

    # Generate random sensor measurements
    u[i] = torch.rand(1, num_sensors)

    # The mapped function equals u the sensors
    # (but we will add a physics-informed loss later w.r.t to all vertices)
    y[i] = x[i]
    v[i] = u[i]

dataset = OperatorDataset(x, u, y, v)

Neural Operator

In this example, we use a DeepONet architecture.

operator = DeepONet(dataset.shapes)

Physics-informed loss

As we want to learn a physically meaningful completion, we define a physics-informed loss function satisfying:

\[- \Delta v = 0\]

y_physics = vertices.clone().requires_grad_(True)[:, ::100] # Use a subset of vertices for the physics-informed loss
mse = torch.nn.MSELoss()

def pde(x, u, y, v):
    # Sensor measurements
    loss = mse(v, u)

    # Physics-informed loss w.r.t. all vertices
    w = lambda y: operator(x, u, y)
    Delta = div(grad(w))

    y_all = y_physics.repeat(u.shape[0], 1, 1).to(y.device)
    loss += (Delta(y_all)**2).mean()

    return loss

loss_fn = PhysicsInformedLoss(pde)

Training

We train the neural operator using the physics-informed loss function.

Trainer(operator, loss_fn=loss_fn).fit(dataset, tol=1e-2, epochs=10000)

Evaluation

Let's evaluate the trained operator on our initial sensor measurements and compute a test error.

x_test = x0.unsqueeze(0)
u_test = u0.unsqueeze(0)

test_loss = mse(operator(x_test, u_test, x_test), u_test)
print(f"loss/test = {test_loss.item():.4e}")

loss/test = 8.9705e-03

Plotting

We can plot the mapped function on the whole mesh by using the cell information from the mesh file and pass it as Triangulation to matplotlib.

tri = Triangulation(vertices[0], vertices[1], mesh.get_cells())

---

Last update: 2024-08-20
Created: 2024-08-20