TaylorGreen

Taylor flow is a two-dimensional analytical solution to the Navier-Stokes equations. The domain is [0, 2] x [0, 2] and the solution is periodic in both x- and y-directions. An interesting feature of the Taylor Green solution is that the pressure and convective terms exactly balance leaving a viscous decay of the initial solution. The analytical velocities and pressure are given at all times as

where is the kinematic viscosity and t is time.

To implement the Taylor Green problem we need to initialize the solution and set up the periodic boundary conditions. The entire problem file TaylorGreen2D.py reads:

from Oasis import *

# Create a mesh here
mesh = RectangleMesh(0, 0, 2, 2, 20, 20)

# Override some problem specific parameters
NS_parameters.update(dict(
    nu = 0.01,         # Kinematic viscosity
    T = 1,             # End time
    dt = 0.001,        # Timestep
    folder = "taylorgreen2D_results",  # folder for storing results
    checkpoint = 1000, # Overwrite solution in Checkpoint folder each checkpoint tstep
    save_step = 10,    # Store solution in new folder each save_step tstep
    max_iter = 1,      # Number of velocity/pressure iterations on each timestep
    use_krylov_solvers = False, 
    velocity_degree = 2,
    pressure_degree = 1,
  )
)

class PeriodicDomain(SubDomain):
    
    def inside(self, x, on_boundary):
        # True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
        return bool((near(x[0], 0) or near(x[1], 0)) and 
              (not ((near(x[0], 0) and near(x[1], 2)) or 
                    (near(x[0], 2) and near(x[1], 0)))) and on_boundary)

    def map(self, x, y):
        if near(x[0], 2) and near(x[1], 2):
            y[0] = x[0] - 2.0
            y[1] = x[1] - 2.0
        elif near(x[0], 2):
            y[0] = x[0] - 2.0
            y[1] = x[1]
        else:
            y[0] = x[0]
            y[1] = x[1] - 2.0

constrained_domain = PeriodicDomain()

initial_fields = dict(
    u0='-sin(pi*x[1])*cos(pi*x[0])*exp(-2.*pi*pi*nu*t)',
    u1='sin(pi*x[0])*cos(pi*x[1])*exp(-2.*pi*pi*nu*t)',
    p='-(cos(2*pi*x[0])+cos(2*pi*x[1]))*exp(-4.*pi*pi*nu*t)/4.')
    
def initialize(q_, q_1, q_2, VV, t, nu, dt, **NS_namespace):
    """Initialize solution. 

    Use t=dt/2 for pressure since pressure is computed in between timesteps.

    """
    for ui in q_:
        deltat = dt/2. if ui is 'p' else 0.
        vv = project(Expression((initial_fields[ui]), t=t+deltat, nu=nu), VV[ui])
        q_[ui].vector()[:] = vv.vector()[:]
        if not ui == 'p':
            q_1[ui].vector()[:] = q_[ui].vector()[:]
            q_2[ui].vector()[:] = q_[ui].vector()[:]

def temporal_hook(q_, t, nu, VV, dt, plot_interval, tstep, **NS_namespace):
    """Function called at end of timestep.    

    Plot solution and compute error by comparing to analytical solution.
    Remember pressure is computed in between timesteps.

    """
    if tstep % plot_interval == 0:
        plot(q_['u0'], title='u')
        plot(q_['u1'], title='v')
        plot(q_['p'], title='p')
    err = {}
    for ui in q_:
        deltat = dt/2. if ui is 'p' else 0.
        vv = project(Expression((initial_fields[ui]), t=t-deltat, nu=nu), VV[ui])
        vv.vector().axpy(-1., q_[ui].vector())
        err[ui] = "{0:2.6f}".format(norm(vv.vector()))
    print "Error at time = ", t, " is ", err

Note that this module contains one mesh and one dictionary called NS_parameters. These are required by the Navier Stokes solver. There are more default parameters found in the NSdefault_hooks.py file. A solution that can be used for visualization will be saved every 10'th timestep (the save_step parameter) to a folder called ´taylorgreen2D_results/Timeseries´. There is also the possibility to store intermediate solutions to a Checkpoint folder. Here we store two previous timesteps of the solution vector and as such the solution can be used to obtain a clean restart (since two previous timesteps are required for the convection term). The solution in the Checkpoint folder is owerwritten continuously, so only the save_step option should be used to create animations etc. The solution is always stored together with the NS_parameters dictionary to allow for clean restarts and to remember solution parameters, like timestep, viscosity etc.

The module also contains a SubDomain named constrained_domain that sets up the two periodic directions. The remaining code is initialization and a function called temporal_hook that is called at the end of each timestep. Here we choose to plot the solution and print out the error by comparing to the analytical solution.

The top of the file NavierStokes.py must be modified to load the correct problem, here by entering the text:

from TaylorGreen2D import *

Running the code in IPython leads to:

In [1]: run NavierStokes
Memory use plain dolfin = 79040
  Inner iterations velocity pressure:
                 error u  error p
    Iter =    1, 1.28e-02 5.56e-02
    Iter =    2, 2.43e-03 1.63e-03
Error at time =  0.001  is  {'p': '0.149662', 'u1': '0.004524', 'u0': '0.004524'}
Error at time =  0.002  is  {'p': '0.150265', 'u1': '0.008806', 'u0': '0.008806'}
Error at time =  0.003  is  {'p': '0.150870', 'u1': '0.012861', 'u0': '0.012861'}
Error at time =  0.004  is  {'p': '0.151598', 'u1': '0.016700', 'u0': '0.016700'}
Error at time =  0.005  is  {'p': '0.152295', 'u1': '0.020335', 'u0': '0.020335'}
Error at time =  0.006  is  {'p': '0.152951', 'u1': '0.023776', 'u0': '0.023776'}
Error at time =  0.007  is  {'p': '0.153571', 'u1': '0.027034', 'u0': '0.027034'}
Error at time =  0.008  is  {'p': '0.154156', 'u1': '0.030119', 'u0': '0.030119'}
Error at time =  0.009  is  {'p': '0.154709', 'u1': '0.033039', 'u0': '0.033039'}
...

The pressure solution plotted in ParaView should look something like:

You may now experiment with higher order elements, shorter timesteps or more velocity/pressure iterations on each timestep.