Fortran file main.f90

This section describes the fortran file main.f90. It focuses on the subroutine my_post_processing that computes the desired outputs. We note that the other subroutines of the file main.f90 should not be modified and are not described here. A template of the file main.f90 is available in the following directory: ($SFEMaNS_DIR)/TEMPLATE.

Remark: If ARPACK is not installed, the following lines of the file main.f90 have to be commented.

1. USE arpack_mhd

2. CALL solver_arpack_mhd(comm_one_d,H_mesh,phi_mesh,&
        inputs%dt,list_mode,mu_H_field)

Apart from these lines, we insist on the fact that only the subroutine my_post_processing should be modified.

This section is splitted into four subsections. First, the structure of the subroutine my_post_processing is described. Then information on functions that compute classic quantities (like L2 norm, divergence, etc.) are given. The third subsection describes how to generate 2D and 3D visualization files for Paraview. Eventually an example is provided.

The subroutine my_post_processing

The subroutine my_post_processing is called after each time iterations. We denote by tn the time after the time iteration n. This subroutine has access to the variables defined at the begining of the file main.f90. Here is a list of the variables that are meaningfull to compute outputs:

1. pp_mesh is the finite element mesh used to approximate the pressure and the level set.
2. vv_mesh is the finite element mesh used to approximate the velocity field.
3. pn is the pressure at time tn.
4. un is the velocity field at time tn.
5. level_set is the level set at time tn (for multiphase computation).
6. density is the density at time tn (for multiphase computation).
7. H_mesh is the finite element mesh used to approximate the magnetic field.
8. phi_mesh is the finite element mesh used to approximate the scalar potential.
9. Hn and Bn are the magnetic fields at time tn. We remind that \(\mu\textbf{H}=\bB\) with \(\mu\) the magnetic permeability.
10. sigma_field is the list of the electrical conductivity of each domains where the magnetic field is approximated.
11. mu_h_field is the list of the magnetic permeability of each domains where the magnetic field is approximated.
12. temp_mesh is the finite element mesh used to approximate the temperature.
13. temperature is the temperature at time tn.
14. temperature_diffusivity_field is the list of the temperature diffusivity of each domains where the temperature is approximated.
15. m_max_c is the number of Fourier modes approximated by each processors.
16. list_mode is the list of the Fourier mode approximated by the processor. We note that m_max_c=SIZE(list_mode).
17. time is the time tn.
18. comm_one_d_ns is the communicator on the Navier-Stokes equations domain.
19. comm_one_d_temp is the communicator on the temperature equation domain.
20. comm_one_d is the communicator on the whole domain. It is equal to the Maxwell equations domain if they are approximated.

The structure of the subroutine my_post_processing is the following:

1. First, the ranks of each processors are defined. It allows to have one processor writing outputs that depends of a specific Fourier mode or region of the finite element mesh. It can be done as follows.

   !===Check ranks
   IF (vv_mesh%me /=0) THEN
      CALL MPI_Comm_rank(comm_one_d_ns(1), rank_ns_S, ierr)
      CALL MPI_Comm_rank(comm_one_d_ns(2), rank_ns_F, ierr)
   ELSE
      rank_ns_s = -1
      rank_ns_f = -1
   END IF
   CALL MPI_Comm_rank(comm_one_d(1), rank_S, ierr)
   CALL MPI_Comm_rank(comm_one_d(2), rank_F, ierr)

   Remarks:

   1. The global rank of each processor, defined earlier in the main.f90, is denoted rank. It is usefull when writing an output that does not depend of a specific Fourier mode or a specific region (like the total kinetic energy).
   2. Comunicators have a dimension equal to two. The first dimension refers to a communication between subdomains of the finite element mesh for a fixed list of Fourier modes. The second dimension refers to a communication between different list of Fourier modes on a fixed subdomain of the finite element mesh.

2. The numerical stability of the computation is checked after each time iteration. It is done as follows.

   !===Check divergence of fields
   IF (inputs%check_numerical_stability) THEN
      IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
         norm = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, un)
      ELSE 
         norm = norm_SF(comm_one_d, 'L2', H_mesh, list_mode, Hn)
      END IF
      IF (norm>1.d2) THEN
         CALL error_petsc('From my_post_processing: numerical unstability')
      END IF
   END IF

   To enable this feature, the following lines are required in the data file:

   ===Check numerical stability (true/false)
   .t.

   It is set to false by default.
3. Verbose and user's outputs are computed every inputs%freq_en time iterations. The parameter inputs%freq_en is set in the data as follows.

   ===Frequency to write restart file
   10

   This part of the subroutine is framed by the following lines code:

   !===Put your postprocessing stuff here
   IF (MOD(it,inputs%freq_en) == 0) THEN
      (verbose and user ouputs)    
   END IF ! end freq_en

4. Generation of 2D and 3D visualization files of the approximated variables. This files can be visualized with the software Paraview. They are generated every inputs%freq_plot time iterations. The parameter inputs%freq_plot is set in the data as follows.

   ===Frequency to create plots
   10

   This part of the subroutine is framed by the following lines code:

   IF (MOD(it,inputs%freq_plot) == 0) THEN
      (3D and 2D plot generation)
   END IF ! end freq_plot

Remark: The first two parts of the subroutine my_post_processing don't need to be modified. Only the section verbose/users ouputs and the section for the visualization files requires modifications depending of the outputs desired.

The verbose ouputs are the average computational time per time iterations, the CFL and the divergence of the velocity and magnetic fields. Here are the required information to have these value computed every inputs%freq_en time iterations.

1. In the verbose/user ouputs section, add the following lines:

   !===Verbose
   CALL write_verbose(rank)

2. Add the following lines in the data file:

   ===Verbose timing? (true/false)
   .t.
   ===Verbose divergence? (true/false)
   .t.
   ===Verbose CFL? (true/false)
   .t.

   Note that computing the divergence does not imply to compute the CFL. These verbose outputs are independent of each other. They are all set to false by default.

Tools to compute user outputs

The subroutine my_post_processing uses the module tn_axi. It gives acces to four functions that can compute various quantity like L2 norm, H1 norm, scalar product of vectors. Here is a description of these functions where we denote by norm a real number.

1. The function dot_product_SF computes the scalar product of two vectors. It is called as follows:

   norm=dot_product_SF(comm_one_d, H_mesh, list_mode, Hn, Bn) 

   where comm_one_d and H_mesh are respectively the communicator and the finite element mesh where the vectors Hn and Bn are defined.
2. The function norm_SF can compute the L2, H1, sH1 norm of a scalar or a vector. It can also computes the L2 norm of the divergence and the curl of a vector. It is called as follows:

   norm = norm_SF(comm_one_d_ns, 'L2', vv_mesh, list_mode, un)

   where comm_one_d_ns and vv_mesh are respectively the communicator and the finite element mesh where un is defined. The characters 'L2' can be switch to H1, sH1, div and curl.
3. The function norm_S can compute the L2, H1, sH1 norm of one Fourier component of a scalar or a vector. It can also computes the L2 norm of the divergence and the curl of one Fourier component of a vector. It is called as follows:

   norm = norm_S(comm_one_d, 'L2', vv_mesh, list_mode(i:i), un(:,:,i:i))

   where comm_one_d_ns and vv_mesh are respectively the communicator and the finite element mesh where un is defined. The characters 'L2' can be switch to H1, sH1, div and curl. The integer i is an integer between 1 and SIZE(list_mode). It represents the label of the Fourier component considered.
4. The function norm_S_L1_zero_mode computes the L1 norm of the Fourier mode zero of a scalar or a vector. It is called as follows:

   norm= norm_S_L1_zero_mode(comm_one_d_ns, pp_mesh, list_mode, pn)

   where comm_one_d_ns and pp_mesh are respectively the communicator and the finite element mesh where pn is defined.

In addition, the file main.f90 contains two subroutines that compute a drag force for problem with solid obstacle and the level set conservation for multiphase problem. The following gives information on these two subroutines.

The subroutine FORCES_AND_MOMENTS computes the drag force of a fluid driven by a vertical velocity field in the presence of a solid obstacle.

1. The presence of the solid obstacle is reported with a penalty function \(\chi\). It is equal to one in the fluid and zero in the solid. The obstacle velocity is denoted \(\bu_\text{obs}\). These functions are both defined in the file condlim.f90. They are only used if the following parameters are set to true in the data file:

   ===Use penalty in NS domain (true/false)?
   .t.
   ===Use nonzero velocity in solids (true/false)?
   .t.

   These parameters are set to false by default. We refer to this section for more information.
2. The drag force can be shown to be equal to \(-\frac{2}{\pi}\int \frac{(1-\chi)(\bu-\bu_\text{obs})\cdot \textbf{e}_z}{dt} \) with \(dt\) being the time step and \(\textbf{e}_z\) being the unit vector in the vertical direction.
3. The inputs are the time, the finite element mesh for the velocity field, the communicator for the Navier-Stokes domain, the list of Fourier mode and the velocity field.
4. The output is written in the file fort.12.
5. This subroutine is called as follows:

   IF (inputs%if_compute_momentum_pseudo_force) THEN
      CALL FORCES_AND_MOMENTS(time,vv_mesh,comm_one_d_ns,list_mode,un)
   END IF

   The parameter inputs%if_compute_momentum_pseudo_force is set to true in the data as follows:

   ===Compute z momentum (true/false)?
   .t.

   The default value is false.

The subroutine compute_level_set_conservation computes the relative error on the level set conservation.

1. For each level set, it computes the term \(\frac{|\int_\Omega (\varphi_{t=tn}-\varphi_{t=0})|}{\int_\Omega \varphi_{t=0}}\).
2. The inputs are the time, the finite element of the level set, the communicator for the Navier-Stokes domain, the list of Fourier mode and the level set.
3. The output is written in the file fort.97.
4. This subroutine is called as follows:

   IF (inputs%if_level_set) THEN
      CALL compute_level_set_conservation(time,pp_mesh,comm_one_d_ns,list_mode,level_set)
   END IF

Tools for visualization with Paraview

The code SFEMaNS can generate 3D and 2D visualization files of real valued scalar function and real valued vector function. These files are generated with the subroutines vtu_3d and make_vtu_file_2D. We note that the 2D plots are generated for each Fourier component of the scalar/vector function.

These subroutines generate files with the extension ".vtu" and ".pvd". The "vtu" files are generated every inputs%freq_plot time iterations. We note that one "vtu" file is generated per processor in meridian section. They contains the informations of the variable to visualize at that time of the computation. One file with the extension ".pvd" is also generated. It contains an overview of all the "vtu" file generated. When opening the "pvd" file with the software Paraview, it allows to have access of the visualization of the variable on the whole domain at different times without loading multiple files.

The file main.f90 provided in the TEMPLATE directory of SFEMaNS generates 3D visualization files for the following variable:

1. the velocity field,
2. the pressure,
3. the level set(s),
4. the density,
5. the temperature,
6. the magnetic field,
7. the scalar potential.

It can also generate 2D files of the above variables when uncommenting the following lines:

!!$    !===VTU 2d======================================================================
!!$    CHARACTER(LEN=3)   :: st_mode
!!$    CHARACTER(LEN=200) :: header
!!$    CHARACTER(LEN=3)   :: name_of_field

and the section between the following lines:

!===Generation of 2D plots for each Fourier mode (uncomment if wanted)
(generation 2D files)
!===End Generation of 2D plots for each Fourier mode (uncomment if wanted)

The subroutines vtu_3d generates 3D visualization files.

1. It is called as follows:

   CALL vtu_3d(var, 'var_mesh', 'File_name', 'Var', what, opt_it=it_plot)

2. The inputs are the following:
   1. var is the variable to plot (like un, pn, Hn, etc.).
   2. var_mesh is the name of the mesh where var is defined (vv_mesh, pp_mesh, temp_mesh, H_mesh or phi_mesh).
   3. File_name is the name of the pvd file generated. The name of the vtu files also start with File_name. However they also present information on the subdomain they represent (_S001, _S002, etc.).
   4. 'Var' is a character of three letters. It is used in Paraview when you want to display the value of the variable var.
   5. what is equal to 'new' or 'old'. If it is 'new', the subroutine generates a pvd file else it only updates the pvd file.
   6. opt_it is an optional integer. It allows to create vtu time every inputs%freq_plot time iterations by completing the name of the vtu files with "_I001", "_I002", etc.
3. The pvd output is named "File_name.pvd". The vtu ouputs have the following format "File_name_S001_I001.vtu".

The subroutine make_vtu_2D generates 2D visualization files.

1. It is called as follows:

   CALL make_vtu_file_2D(comm_one_d_var(1), var_mesh, header, var(:,:,i), name_of_field, what, opt_it=it_plot)

2. The inputs are the following:
   1. comm_one_d_var(1) is the comunicator associated to the variable var. Only the first dimension is given. It represents communication between subsections of the finite element mesh (the Fourier component is fixed).
   2. var_mesh is the name of the mesh where var is defined (vv_mesh, pp_mesh, temp_mesh, H_mesh or phi_mesh).
   3. header is the name of the pvd file. The name of the vtu files also start with header. However they also present information on the subdomain they represent (_S001, _S002, etc.).
   4. var(:,:,i) is Fourier component of the variable to plot. We note that the integer i is only the label of the Fourier component, the Fourier mode is equal to list_mode(i).
   5. name_of_field is a character of three letters. It is used in Paraview when you want to display the value of the variable var.
   6. what is equal to 'new' or 'old'. If it is new, the subroutine generates a pvd file else it only updates the pvd file.
   7. opt_it is an optional integer. It allows to create vtu time every inputs%freq_plot time iterations by completing the name of the vtu files with "_I001", "_I002", etc.
3. The pvd output is named "header.pvd". The vtu ouputs have the following format "header_S001_I001.vtu".

Example

We give a description of the subroutine my_post_processing of the template file main.f90. This file can be found in the following directory: ($SFEMaNS_DIR)/TEMPLATE. We refer to the section Examples on physical problems for more examples.

1. Modules are added at the begining of the subroutine my_post_processing.

   USE sub_plot
   USE chaine_caractere
   USE tn_axi
   USE boundary
   USE sft_parallele
   USE verbose

2. The declaration of local variables is done as follows.

       IMPLICIT NONE
       INTEGER,                             INTENT(IN) :: it
       REAL(KIND=8)                                    :: err, norm
       INTEGER                                         :: i, it_plot
       CHARACTER(LEN=3)                                :: what
       INTEGER                                         :: rank_S, rank_F 
       INTEGER                                         :: rank_ns_S, rank_ns_F
       REAL(KIND=8), DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: level_1
   !!$    !===VTU 2d======================================================================
   !!$    CHARACTER(LEN=3)   :: st_mode
   !!$    CHARACTER(LEN=200) :: header
   !!$    CHARACTER(LEN=3)   :: name_of_field

3. The ranks of each processors are defined.

   !===Check ranks
   IF (vv_mesh%me /=0) THEN
      CALL MPI_Comm_rank(comm_one_d_ns(1), rank_ns_S, ierr)
      CALL MPI_Comm_rank(comm_one_d_ns(2), rank_ns_F, ierr)
   ELSE
      rank_ns_s = -1
      rank_ns_f = -1
   END IF
   CALL MPI_Comm_rank(comm_one_d(1), rank_S, ierr)
   CALL MPI_Comm_rank(comm_one_d(2), rank_F, ierr)

4. The numerical stability of the computation is checked after each time iteration.

   !===Check divergence of fields
   IF (inputs%check_numerical_stability) THEN
      IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
         norm = norm_SF(comm_one_d_ns, 'L2', vv_mesh, list_mode, un)
      ELSE 
         norm = norm_SF(comm_one_d, 'L2', H_mesh, list_mode, Hn)
      END IF
      IF (norm>1.d2) THEN
         CALL error_petsc('From my_post_processing: numerical unstability')
      END IF
   END IF

5. The verbose and user's outputs are computed every inputs%freq_en time iterations.

   !===Put your postprocessing stuff here
   IF (MOD(it,inputs%freq_en) == 0) THEN

   1. The verbose outputs are computed as follows.

      !===Verbose
      CALL write_verbose(rank)

   2. The user ouputs that involve the Navier-Stokes equations and the temperature equation are computed if these equations are approximated. It is done by checking the type of problem approximated.

      IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN

      1. If the parameter inputs%if_compute_momentum_pseudo_force is set to true is the data, the drag force is computed. It is written in the file fort.12.

         IF (inputs%if_compute_momentum_pseudo_force) THEN
            !===Compute the term -2/pi*integral((1-chi)*(u-u_solid)/dt)
            !===chi is the penalty function, u_solid the velocity of the solid, dt the time step
            !==Output written in the file fort.12
            CALL FORCES_AND_MOMENTS(time,vv_mesh,comm_one_d_ns,list_mode,un)
         END IF

      2. The divergence of the velocity field is computed as follows. It is written is the file fort.31.

         err = norm_SF(comm_one_d, 'div', vv_mesh, list_mode, un)
         norm = norm_SF(comm_one_d, 'H1', vv_mesh, list_mode, un)
         IF (rank == 0) THEN
            !===Divergence of velocity field
            WRITE(31,*) time, err/norm
         END IF

      3. The L2 norm of each Fourier component m of the velocity field is written is the file fort.xxx with xxx=100+m. It uses the function norm_S as follows.

         DO i=1,SIZE(list_mode)
            norm = norm_S(comm_one_d, 'L2', vv_mesh, list_mode(i:i), un(:,:,i:i))
            IF (rank_ns_S == 0) THEN
               !===L2 norm of Fourier mode list_mode(i) of velocity field un
               WRITE(100+list_mode(i),*) time, norm 
            END IF
         END DO

      4. The L2 norm of the velocity is computed as follows. It is written in the file fort.98.

         err = norm_SF(comm_one_d, 'L2', vv_mesh, list_mode, un)
         norm = norm_SF(comm_one_d, 'sH1', vv_mesh, list_mode, un)
         IF (rank == 0) THEN
            WRITE(98,*) 'norm L2 of velocity', time, err
            WRITE(*,*) 'norm L2 of velocity', time, err
            WRITE(*,*) 'semi norm H1 of velocity', time, norm
         END IF

      5. The L2 norm of the pressure is computed as follows.

         err = norm_SF(comm_one_d, 'L2', pp_mesh, list_mode, pn)
         IF (rank == 0) THEN
            WRITE(*,*) 'norm L2 of pressure', time, err
         END IF

      6. If inputs%if_level_set is set to true (multiphase problem), the conservation of the level set is computed as follows. It is written in the file fort.97.

         IF (inputs%if_level_set) THEN
            !===Compute the term integral(level_set-level_set_t=0)/integral(level_set_t=0)
            !===Output written in file fort.97
            CALL compute_level_set_conservation(time,pp_mesh,comm_one_d_ns,list_mode,level_set)
         END IF

      7. If inputs%if_temperature is set to true, the L2 norm of the temperature is computed as follows.

         IF (inputs%if_temperature) THEN
            err = norm_SF(comm_one_d, 'L2', temp_mesh, list_mode, temperature)
            IF (rank == 0) THEN
               WRITE(*,*) 'norm L2 of temperature', time, err
            END IF
         END IF

   3. End of the user outputs involving the Navier-Stokes equations and the temperature equation.

      END IF ! end nst or mhd or fhd

   4. The user outputs involving the Maxwell equations are computed if these equations are approximated. It is done by checking the type of problem approximated.

      IF (inputs%type_pb/='nst') THEN

      1. The L2 norm of the magnetic field \(\bH\) is computed as follows. It is written in the file fort.41.

         err = norm_SF(comm_one_d, 'L2', H_mesh, list_mode, Hn)
         IF (rank == 0) THEN
            !===L2 norm of magnetic field
            WRITE(41,*) time, err
         END IF

      2. The divergence and the L2 norm of the magnetic field \(\bB\) are computed as follows. They are written in the files fort.51 and fort.52.

         err = norm_SF(comm_one_d, 'div', H_mesh, list_mode, Bn)
         norm = norm_SF(comm_one_d, 'L2', H_mesh, list_mode, Bn)
         IF (rank == 0) THEN
            !===L2 norm of div(Bn)
            WRITE(51,*) time, err, err/norm
            WRITE(52,*) time, err, norm
            WRITE(*,*) 'norm L2 of magnetic field', time, norm
         END IF

      3. The L2 norm of each Fourier component m of the magnetic field \(\bH\) is written is the file fort.xxx with xxx=200+m. It uses the function norm_S as follows.

         DO i=1,SIZE(list_mode)
            norm = norm_S(comm_one_d, 'L2', H_mesh, list_mode(i:i), Hn(:,:,i:i))
            IF (rank_S == 0) THEN
               !===L2 norm of Fourier mode list_mode(i) of magnetic field Hn
               WRITE(200+list_mode(i),*) time, norm
            END IF
         END DO

   5. End of the user ouputs involving the Maxwell equations.

      END IF ! end /=nst

6. End of the verbose and user ouputs.

   END IF ! end freq_en

7. The visualization files are generated every inputs%freq_plot time iterations.

   IF (MOD(it,inputs%freq_plot) == 0) THEN

   1. The inputs what and it_plot of the subroutines vtu_3d and make_vtu_file_2D are defined as follows.

      !===Plot whatever you want here
      IF (it==inputs%freq_plot) THEN
         what = 'new'
      ELSE
         what = 'old'
      END IF
      it_plot = it/inputs%freq_plot

   2. The generation of the 3D files is done with the subroutine vtu_3d as follows.

      !===Generation of 3D plots
      IF (inputs%if_level_set) THEN
         level_1=level_set(1,:,:,:)
         CALL vtu_3d(density, 'vv_mesh', 'Density', 'density', what, opt_it=it_plot)
         CALL vtu_3d(level_1, 'pp_mesh', 'Level_1', 'level_1', what, opt_it=it_plot)
         IF (SIZE(level_set,1).GE.2) THEN
            level_1=level_set(2,:,:,:)
            CALL vtu_3d(level_1, 'pp_mesh', 'Level_2', 'level_2', what, opt_it=it_plot)
         END IF
      END IF
      IF (inputs%type_pb/='mxw' .AND. inputs%type_pb/='mxx') THEN
         CALL vtu_3d(un, 'vv_mesh', 'Velocity', 'vel', what, opt_it=it_plot)
         CALL vtu_3d(pn, 'pp_mesh', 'Pressure', 'pre', what, opt_it=it_plot)
      END IF
      IF (inputs%if_temperature) THEN
         CALL vtu_3d(temperature, 'temp_mesh', 'Temperature', 'temp', what, opt_it=it_plot)
      END IF
      IF (inputs%type_pb/='nst') THEN
         CALL vtu_3d(Hn, 'H_mesh', 'MagField', 'mag', what, opt_it=it_plot)
         IF (inputs%nb_dom_phi>0) THEN
            CALL vtu_3d(phin, 'phi_mesh', 'ScalPot', 'phi', what, opt_it=it_plot)
         END IF
      END IF
      !==End generation of 3D plots

   3. The generation of the 2D files is commented. It uses the subroutine make_vtu_file_2D.

             !===Generation of 2D plots for each Fourier mode (uncomment if wanted)
             !===Proceed as follows to make 2D plots in the Fourier space (using Hn for instance)
             !===what = 'new' if first image, what= 'old' for later images  (CHARACTER(LEN=3)   :: what)
             !===WRITE(st_mode,'(I3)') list_mode(i)                         (CHARACTER(LEN=3)   :: st_mode)
             !===header = 'Hn_'//'_mode_'//trim(adjustl(st_mode))           (CHARACTER(LEN=200) :: header)
             !===name_of_field = 'Hn' (for instance)                        (CHARACTER(LEN=3)   :: name_of_field)
             !===CALL make_vtu_file_2D(comm_one_(1), H_mesh, header, Hn(:,:,i), name_of_field, what, opt_it=1)
      !!$       IF (inputs%if_level_set) THEN
      !!$          !===2D plots for each mode of the first level set
      !!$          DO i = 1, m_max_c
      !!$             WRITE(st_mode,'(I3)') list_mode(i)
      !!$             header = 'Ln_'//'mode_'//trim(adjustl(st_mode)) 
      !!$             name_of_field = 'Ln'
      !!$             CALL make_vtu_file_2D(comm_one_d_ns(1), pp_mesh, header, level_1(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$             header = 'Dn_'//'mode_'//trim(adjustl(st_mode))
      !!$             name_of_field = 'Dn'
      !!$             CALL make_vtu_file_2D(comm_one_d_ns(1), vv_mesh, header, density(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$          END DO
      !!$       END IF
      !!$       IF (inputs%type_pb/='mxw' .AND. inputs%type_pb/='mxx') THEN
      !!$          !===2D plots for each mode of the velocity field and the pressure
      !!$          DO i = 1, m_max_c
      !!$             WRITE(st_mode,'(I3)') list_mode(i)
      !!$             header = 'Vn_'//'mode_'//trim(adjustl(st_mode))
      !!$             name_of_field = 'Vn'
      !!$             CALL make_vtu_file_2D(comm_one_d_ns(1), vv_mesh, header, un(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$             WRITE(st_mode,'(I3)') list_mode(i)
      !!$             header = 'Pn_'//'mode_'//trim(adjustl(st_mode))
      !!$             name_of_field = 'Pn'
      !!$             CALL make_vtu_file_2D(comm_one_d_ns(1), pp_mesh, header, pn(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$          END DO
      !!$       END IF
      !!$       IF (inputs%if_temperature) THEN
      !!$          !===2D plots for each mode of the temperature
      !!$          DO i = 1, m_max_c
      !!$             WRITE(st_mode,'(I3)') list_mode(i)
      !!$             header = 'Tn_'//'_mode_'//trim(adjustl(st_mode))
      !!$             name_of_field = 'Tn'
      !!$             CALL make_vtu_file_2D(comm_one_d(1), temp_mesh, header, temperature(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$          END DO
      !!$       END IF
      !!$       IF (inputs%type_pb/='nst') THEN
      !!$          !===2D plots for each mode of the magnetic field and the scalar potential
      !!$          DO i = 1, m_max_c
      !!$             WRITE(st_mode,'(I3)') list_mode(i)
      !!$             header = 'Hn_'//'_mode_'//trim(adjustl(st_mode))
      !!$             name_of_field = 'Hn'
      !!$             CALL make_vtu_file_2D(comm_one_d(1), H_mesh, header, Hn(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$             IF (inputs%nb_dom_phi>0) THEN
      !!$                WRITE(st_mode,'(I3)') list_mode(i)
      !!$                header = 'Phin_'//'_mode_'//trim(adjustl(st_mode))
      !!$                name_of_field = 'Phin'
      !!$                CALL make_vtu_file_2D(comm_one_d(1), phi_mesh, header, phin(:,:,i), name_of_field, what, opt_it=it_plot)
      !!$             END IF
      !!$          END DO
      !!$       END IF
             !===End Generation of 2D plots for each Fourier mode (uncomment if wanted)

8. End of the section that generates visualization files.

   END IF ! end freq_plot

Remark: If more than 100 Fourier modes are used, a conflict arises between different energy files. Indeed, the L2 norm of the Fourier component i of the velocity is written in the files fort.xxx with xxx=100+i. The same is done for the magnetic field with the files fort.yyy with yyy=200+i. This problem can be overcome by switching 200 to a larger number. One can also defines proper name for the energy outputs file.