Model composition

Table of Contents

• Problem definition
• Input file
• Lauching the projectiles

Problem definition

Recall that the complete model for simulating the projectile trajectory is

\begin{align} \dot{\boldsymbol{x}} & = \boldsymbol{v}, \label{1} \\ \dot{\boldsymbol{v}} & = \boldsymbol{a} = \boldsymbol{g} - \mu \boldsymbol{v}, \label{2} \\ \mathbf{r} = & = \begin{Bmatrix} \tilde{\boldsymbol{x}} - \boldsymbol{x}_n - \left(t - t_n\right) \dot{\boldsymbol{x}} \\ \tilde{\boldsymbol{v}} - \boldsymbol{v}_n - \left(t - t_n\right) \dot{\boldsymbol{v}} \\ \end{Bmatrix}, \label{3} \\ \begin{Bmatrix} \boldsymbol{x} \\ \boldsymbol{v} \end{Bmatrix} & = \mathop{\mathrm{root}}\limits_{\tilde{\boldsymbol{x}}, \tilde{\boldsymbol{v}}} \left( \mathbf{r} \right), \label{4} \end{align}

subject to appropriate initial conditions \(\boldsymbol{x}_0\) and \(\boldsymbol{v}_0\)

Among these equations:

• \(\eqref{1}\) and \(\eqref{3}\) can be defined using VecBackwardEulerTimeIntegration.
• \(\eqref{2}\) is the custom model ProjectileAcceleration which we have implemented in previous tutorials.
• \(\eqref{4}\) is the ImplicitUpdate.

This tutorial demonstrates that our custom model ProjectileAcceleration can be composed with other existing, predefined NEML2 models.

Input file

The following input file composes the constitutive model for a single-step update for 5 projectiles each with a different dynamic viscosity, i.e., the shape of "mu" is \((5;)\).

[Tensors]
  [g]
    type = Vec
    values = '0 -9.81 0'
  []
  [mu]
    type = Scalar
    values = '0.01 0.05 0.1 0.5 1'
    batch_shape = (5)
  []
[]
 
[Models]
  [eq2]
    type = ProjectileAcceleration
    velocity = 'state/v'
    acceleration = 'state/a'
    gravitational_acceleration = 'g'
    dynamic_viscosity = 'mu'
  []
  [eq3a]
    type = VecBackwardEulerTimeIntegration
    variable = 'state/x'
    rate = 'state/v'
  []
  [eq3b]
    type = VecBackwardEulerTimeIntegration
    variable = 'state/v'
    rate = 'state/a'
  []
  [eq3]
    type = ComposedModel
    models = 'eq3a eq3b'
  []
  [system]
    type = ComposedModel
    models = 'eq2 eq3'
  []
  [eq4]
    type = ImplicitUpdate
    implicit_model = 'system'
    solver = 'newton'
  []
[]
 
[Solvers]
  [newton]
    type = Newton
    rel_tol = 1e-08
    abs_tol = 1e-10
    max_its = 50
  []
[]

Lauching the projectiles

To obtain the entire trajectory, the constitutive model need to be recursively integrated. As discussed in a previous tutorial, a transient driver should be used to perform the recursive constitutive update. Conveniently, in this example, since we are dealing with an autonomous system (i.e., no external "forces" are needed to drive the constitutive update), the vanilla TransientDriver can be used:

[Drivers]
  [driver]
    type = TransientDriver
    model = 'eq4'
    prescribed_time = 'times'
    ic_Vec_names = 'state/x state/v'
    ic_Vec_values = 'x0 v0'
    save_as = 'result.pt'
    show_input_axis = true
    show_output_axis = true
  []
[]
 
[Tensors]
  [end_time]
    type = Scalar
    values = '1'
    batch_shape = (1,1)
  []
  [times]
    type = LinspaceScalar
    start = 0
    end = 'end_time'
    nstep = 100
  []
  [x0]
    type = Vec
    values = '0 0 0'
  []
  [v0]
    type = Vec
    values = "10 5 0
              8 6 0"
    batch_shape = (2,1)
  []
[]

The five projectiles are launched from the same position (the origin) but with two different lauching velocities. Note how broadcasting is used to simultaneously simulate the trajectories of all 10 combinations.

#include "neml2/drivers/Driver.h"
 
int
main()
{
  using namespace neml2;
  set_default_dtype(kFloat64);
  auto factory = load_input("input.i");
  auto driver = factory->get_driver("driver");
  driver->run();
}

Output

eq4's input axis:
  0: forces/t     [(0, 1)]
  1: old_forces/t [(1, 2)]
  2: old_state/v  [(2, 5)]
  3: old_state/x  [(5, 8)]
  4: state/v      [(8, 11)]
  5: state/x      [(11, 14)]
eq4's output axis:
  0: state/v [(0, 3)]
  1: state/x [(3, 6)]

The following Python script plots the trajectories loaded from "result.pt" written by the driver.

import torch
from matplotlib import pyplot as plt
 
# Load the result
res = torch.jit.load("result.pt")
O = dict(res.output.named_buffers())
 
# Plot
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
nstep = 100
nvel = 2 # two launching velocities
nproj = 5 # five projectiles (with different dynamic viscosity)
for i in range(nvel):
    for j in range(nproj):
        x = [O["{}.state/x".format(n)][i, j, 0].item() for n in range(1, nstep)]
        y = [O["{}.state/x".format(n)][i, j, 1].item() for n in range(1, nstep)]
        ax[i].plot(x, y, "--")
    ax[i].set_xlabel("x")
    ax[i].set_ylabel("y")
    ax[i].grid()
fig.tight_layout()
fig.savefig("trajectories.svg")

Projectile trajectories

Previous	Next
The forward operator	Model visualization