The subroutine init_velocity_pressure

It is used to initialize the velocity field, the pressure and the pressure increment.

Inputs and outputs

The inputs of this function are the following:

1. mesh_f is the finite element mesh used to approximate the velocity field.
2. mesh_c is the finite element mesh used to approximate the pressure.
3. dt is the time step.
4. list_mode is a list of integers which contains the Fourier modes approximated.

As the mesh can be subdivised, we note that mesh_f and mesh_c depend of the processor considered when doing parallel computing. Same goes for the list of Fourier mode list_mode.

The outputs of this function are the following:

1. time is the time when the computations starts.
2. un_m1 is the velocity field at the time time-dt.
3. un is the velocity field at the time time.
4. pn_m1 is the pressure at the time time-dt.
5. pn is the pressure at the time time.
6. phin_m1 is the increment pressure at the time time-dt.
7. phin is the increment pressure at the time time.

Remarks:

1. For a velocity field with a zero divergence, the pressure increments satisfy the relations: phin=pn-pn_m1 and phin_m1=pn_m1-pn_m2.
2. The velocity field is a vector so its format is a real valued tabular of three columns with dimension (mesh_f%np,6,SIZE(list_mode)). We remind that mesh_f%np is the number of nodes of the finite element mesh mesh_f.
3. The pressure and increment pressure are scalar so their format is a real valued tabular of three columns of dimension (mesh_c%np,2,SIZE(list_mode)). We remind that mesh_c%np is the number of nodes of the finite element mesh mesh_p.

Exemple

Here is an exemple where we use the function vv_exact and pp_exact to initialize the velocity field and the pressire. Thus, the initial conditions satisfy the boundary conditions.

time = 0.d0
 DO i= 1, SIZE(list_mode)
    mode = list_mode(i) 
    DO j = 1, 6 
       !===velocity
       un_m1(:,j,i) = vv_exact(j,mesh_f%rr,mode,time-dt)  
       un   (:,j,i) = vv_exact(j,mesh_f%rr,mode,time)
    END DO
    DO j = 1, 2
       !===pressure
       pn_m2(:)       = pp_exact(j,mesh_c%rr,mode,time-2*dt)
       pn_m1  (:,j,i) = pp_exact(j,mesh_c%rr,mode,time-dt)
       pn     (:,j,i) = pp_exact(j,mesh_c%rr,mode,time)
       phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
       phin   (:,j,i) = Pn   (:,j,i) - pn_m1(:,j,i)
    ENDDO
 ENDDO
 RETURN

We note that the integers mode, i and j have to be declared. Same for the real valued tabular pn_m2 whose dimension is mesh_c%np. It is done by adding the two following lines in the declaration of the function init_velocity_pressure:

INTEGER                                    :: mode, i, j 
REAL(KIND=8), DIMENSION(mesh_c%np)         :: pn_m2

We refer to the sections Examples with manufactured solutions and Examples on physical problems for more examples.