These Notes Describe Code Modifications Done During the Lectures, Concerning

Notes aboutcodemodificationsduringtheLectures

These notes describe code modifications done during the lectures, concerning:

a)Inclusion of Central Differences Advection,

b)Implicitcalculation.

The initial code was using upwind for advection and an explicit calculation. In the case the coefficients of the equation are:

Central Differences for Advection

This chapter shows how to compute central differences, both in explicit and implicit options:

In case of Central-differences discretized explicit equation is:

Central differences explicit

If calculation is explicit one gets:

That can be written in the canonic for is:

With the following coefficients:

In the code, these coefficients are stored in arrays coef.left(i), coef.center(i) and coef.right(i).

In case of Explicit calculation those coefficients are:

Central differences implicit

If Calculation is implicit one would get:

Moving the unknowns to le left member one would get

Or, in the canonic form:

In the code, these coefficients are stored in arrays coef.left(i), coef.center(i) and coef.right(i) and coef.Ti(I).

In case of implicit calculation those coefficients are:

Implicit calculation

In case of implicit calculation, one has 3 unknowns in each equation and as many equations as grid points and consequently it is necessary to invert a matrix. The figure below shows how to organize the matrix. On columns one sets the unknowns, meaning that we need as many columns as grid points. On the lines one has to write the equations, one for each grid point. Like this we get a square matrix. In each line one has 3 coefficients, one multiplying by the local unknown (the main diagonal, center.coef(I)), one multiplying by the unknown on the left (lower diagonal, left.coef(I)) and another multiplying by the value on the right (upper diagonal, right.coef(I)). The independent term is in case of the transport the old value property value (ti.coef(I))

One can easily see that in first equation the (grid point 1) there is no room for the left.coef(1). In fact, this coefficient will multiply by the C0 that must be known, i.e. it is a boundary condition. Being known, the multiplication can be done and the value can be moved to the independent term:

And identically for the upper diagonal for the last point (NC). In that case one gets:

Being a tridiagonal matrix, it is very easy to invert.

The inversion can be done using the Thomas algorithm or the Gauss elimination using a “forward elimination->backward substitution”. Defining the multiplier F as:

And adding line (i-1) multiplied by F to the line (i) one can eliminate the lower diagonal if we start with the second line. Doing so one can compute the value of C at the NC cell directly. After that one can perform a backward substitution to compute the other values of C.