Moving Time Window: Taylor Green Vortex Decay

Introduction

Some of the previous tutorials have shown solving the transient problems with continuous time approach. This tutorial presents a moving time window approach to solve a complex transient Navier-Stokes problem. In this tutorial, you will learn:

  1. How to solve sequences of problems/domains in Modulus.

  2. How to set up periodic boundary conditions.

Note

This tutorial assumes you have completed tutorial on 1D Wave Equation.

Problem Description

As mentioned in tutorial 1D Wave Equation, solving transient simulations with only the continuous time method can be difficult for long time durations. This tutorial shows how this can be overcome by using a moving time window. The example problem is the 3D Taylor-Green vortex decay at Reynolds number 500. The Taylor-Green vortex problem is often used as a benchmark to compare solvers and in this case we will generate validation data with a spectral solver. The domain is a cube of length \(2\pi\) with periodic boundary conditions on all sides. The initial conditions you will use are,

(178)\[\begin{split}\begin{aligned} u(x,y,z,0) &= sin(x)cos(y)cos(z) \\ v(x,y,z,0) &= -cos(x)sin(y)cos(z) \\ w(x,y,z,0) &= 0 \\ p(x,y,z,0) &= \frac{1}{16} (cos(2x) + cos(2y))(cos(2z) + 2) \\\end{aligned}\end{split}\]

This problem shows solving the time dependent incompressible Navier-Stokes equations with a density \(1\) and viscosity \(0.002\). Note that because there are periodic boundaries on all sides you do not need boundary conditions Fig. 90.

Taylor-Green vortex initial conditions.

Fig. 90 Taylor-Green vortex initial conditions.

The moving time window approach works by iteratively solving for small time windows to progress the simulation forward. The time windows use the previous window as new initial conditions. The continuous time method is used for solving inside a particular window. A figure of this method can be found here Fig. 91 for a hypothetical 1D problem. Learning rate decay is restarted after each time window. A similar approach can be found here 1.

Moving Time Window Method

Fig. 91 Moving Time Window Method

Case Setup

The case setup for this problem is similar to many of the previous tutorials except for two key differences. This example shows how to set up a sequence of domains to iteratively solve for.

Note

The python script for this problem can be found at examples/taylor_green/.

Sequence of Train Domains

First, construct your geometry similar to the previous problems.

Also, define values for how large our time window will be. In this case solve to \(1\) unit of time and then construct all needed nodes.

    # time window parameters
    time_window_size = 1.0
    t_symbol = Symbol("t")
    time_range = {t_symbol: (0, time_window_size)}
    nr_time_windows = 10

    # make navier stokes equations
    ns = NavierStokes(nu=0.002, rho=1.0, dim=3, time=True)

    # define sympy variables to parametrize domain curves
    x, y, z = Symbol("x"), Symbol("y"), Symbol("z")

    # make geometry for problem
    channel_length = (0.0, 2 * np.pi)
    channel_width = (0.0, 2 * np.pi)
    channel_height = (0.0, 2 * np.pi)
    box_bounds = {x: channel_length, y: channel_width, z: channel_height}

    # define geometry
    rec = Box(
        (channel_length[0], channel_width[0], channel_height[0]),
        (channel_length[1], channel_width[1], channel_height[1]),
    )

The architecture can be created as below.

    # make network for current step and previous step
    flow_net = FullyConnectedArch(
        input_keys=[Key("x"), Key("y"), Key("z"), Key("t")],
        output_keys=[Key("u"), Key("v"), Key("w"), Key("p")],
        periodicity={"x": channel_length, "y": channel_width, "z": channel_height},
        layer_size=256,
    )
    time_window_net = MovingTimeWindowArch(flow_net, time_window_size)

Note that the periodicity is set when creating out network architecture. This will force the network by construction to give periodic solutions on the given bounds. Now two domains are defined, one for the initial conditions and one to be used for all future time windows.

    # make initial condition domain
    ic_domain = Domain("initial_conditions")

    # make moving window domain
    window_domain = Domain("window")

    # make initial condition
    ic = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=rec,
        outvar={
            "u": sin(x) * cos(y) * cos(z),
            "v": -cos(x) * sin(y) * cos(z),
            "w": 0,
            "p": 1.0 / 16 * (cos(2 * x) + cos(2 * y)) * (cos(2 * z) + 2),
        },
        batch_size=cfg.batch_size.initial_condition,
        bounds=box_bounds,
        lambda_weighting={"u": 100, "v": 100, "w": 100, "p": 100},
        param_ranges={t_symbol: 0},
    )
    ic_domain.add_constraint(ic, name="ic")

    # make constraint for matching previous windows initial condition
    ic = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=rec,
        outvar={"u_prev_step_diff": 0, "v_prev_step_diff": 0, "w_prev_step_diff": 0},
        batch_size=cfg.batch_size.interior,
        bounds=box_bounds,
        lambda_weighting={
            "u_prev_step_diff": 100,
            "v_prev_step_diff": 100,
            "w_prev_step_diff": 100,
        },
        param_ranges={t_symbol: 0},
    )
    window_domain.add_constraint(ic, name="ic")

    # make interior constraint
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=rec,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
        bounds=box_bounds,
        batch_size=4094,
        param_ranges=time_range,
    )
    ic_domain.add_constraint(interior, name="interior")
    window_domain.add_constraint(interior, name="interior")

    # add inference data for time slices
    for i, specific_time in enumerate(np.linspace(0, time_window_size, 10)):
        vtk_obj = VTKUniformGrid(
            bounds=[(0, 2 * np.pi), (0, 2 * np.pi), (0, 2 * np.pi)],
            npoints=[128, 128, 128],
            export_map={"u": ["u", "v", "w"], "p": ["p"]},
        )
        grid_inference = PointVTKInferencer(
            vtk_obj=vtk_obj,
            nodes=nodes,
            input_vtk_map={"x": "x", "y": "y", "z": "z"},
            output_names=["u", "v", "w", "p"],
            requires_grad=False,
            invar={"t": np.full([128 ** 3, 1], specific_time)},
            batch_size=100000,
        )
        ic_domain.add_inferencer(grid_inference, name="time_slice_" + str(i).zfill(4))
        window_domain.add_inferencer(
            grid_inference, name="time_slice_" + str(i).zfill(4)
        )

In the moving time window domain there is a new initial condition that will come from the previous time window. Now that the domain files have been constructed for the initial conditions and future time windows You can put it all together and make your sequential solver.

Sequence Solver

    # make solver
    slv = SequentialSolver(
        cfg,
        [(1, ic_domain), (nr_time_windows, window_domain)],
        custom_update_operation=time_window_net.move_window,
    )

    # start solver
    slv.solve()

Instead of using the normal Solver class use the SequentialSolver class. This class expects a list where each element is a tuple of the number of times solving as well as the given domain. In this case, solve the initial condition domain once and then solve the time window domain multiple times. A custom_update_operation is provided which is called after solving each iteration and in this case is used to update network parameters.

The iterative domains are solved here in a very general way. In subsequent chapters this structure is used to implement a complex iterative algorithm to solve conjugate heat transfer problems.

Results and Post-processing

After solving this problem you can visualize the results. Look in the network directory to see,

current_step.txt initial_conditions window_0000 window_0001...

If you look in any of these directories you see the typical files storing network checkpoint and results. Plotting at time 15 gives the snapshot Fig. 90. To validate these results we have conducted the same simulation using a spectral solver. Comparing the point-wise error of transient simulations can be misleading as this is a chaotic dynamical system and any small errors will quickly cause large differences. Instead, look at the average turbulent kinetic energy decay (TKE) of the simulation. A figure of this can be found here below.

Taylor-Green vortex at time :math:`15.0`.

Fig. 92 Taylor-Green vortex at time \(15.0\).

Taylor-Green Trubulent kinetic energy decay.

Fig. 93 Taylor-Green Trubulent kinetic energy decay.

References

1

Colby L Wight and Jia Zhao. Solving allen-cahn and cahn-hilliard equations using the adaptive physics informed neural networks. arXiv preprint arXiv:2007.04542, 2020.