sub_temperature.f90

Go to the documentation of this file.

!
!Authors Jean-Luc Guermond, Raphael Laguerre, Caroline Nore, Copyrights 2005
!Revised for PETSC, Jean-Luc Guermond, Franky Luddens, January 2011
!
MODULE subroutine_temperature
  USE my_util
  USE boundary

  PUBLIC :: three_level_temperature
  PRIVATE
CONTAINS

  SUBROUTINE three_level_temperature(comm_one_d,time, temp_1_LA, dt, list_mode, &
       temp_mesh, tempn_m1, tempn, chmp_vit, vol_heat_capacity, &
       temp_diffusivity, my_par_cc, temp_list_dirichlet_sides, temp_per)
    !============================== 
    USE def_type_mesh
    USE fem_m_axi
    USE fem_rhs_axi
    USE fem_tn_axi
    USE dir_nodes_petsc
    USE periodic
    USE st_matrix
    USE solve_petsc
    USE dyn_line
    USE boundary
    USE chaine_caractere
    USE sub_plot
    USE st_matrix
    USE sft_parallele
    IMPLICIT NONE
    REAL(KIND=8)                                        :: time, dt
    INTEGER,      DIMENSION(:),     INTENT(IN)          :: list_mode 
    TYPE(mesh_type),                INTENT(IN)          :: temp_mesh
    type(petsc_csr_la)                                  :: temp_1_la
    TYPE(periodic_type),            INTENT(IN)          :: temp_per
    TYPE(solver_param),             INTENT(IN)          :: my_par_cc
    REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT)       :: tempn_m1, tempn
    INTEGER,      DIMENSION(:),     INTENT(IN)          :: temp_list_dirichlet_sides 
    REAL(KIND=8), DIMENSION(:),     INTENT(IN)          :: vol_heat_capacity, temp_diffusivity
    REAL(KIND=8), DIMENSION(:,:,:),          INTENT(IN) :: chmp_vit
    LOGICAL,                                       SAVE :: once = .true.
    INTEGER,                                       SAVE :: m_max_c
    INTEGER,     DIMENSION(:),   POINTER,          SAVE :: temp_js_d ! Dirichlet nodes
    INTEGER,                                       SAVE :: my_petscworld_rank
    REAL(KIND=8),                                  SAVE :: mass0, hmoy
    !----------FIN SAVE--------------------------------------------------------------------

    !----------Declaration sans save-------------------------------------------------------
    INTEGER,          POINTER, DIMENSION(:)  :: temp_1_ifrom
    INTEGER                                  :: i, m, n, l, index
    INTEGER                                  :: code, mode
    !allocations des variables locales
    REAL(KIND=8), DIMENSION(temp_mesh%np)                      :: ff
    REAL(KIND=8), DIMENSION(temp_mesh%np, 2)                   :: tempn_p1
    REAL(KIND=8), DIMENSION(temp_mesh%gauss%l_G*temp_mesh%me,2, SIZE(list_mode)) :: ff_conv
    REAL(KIND=8)   ::tps, tps_tot, tps_cumul
    REAL(KIND=8) :: one, zero, three
    DATA zero, one, three/0.d0,1.d0,3.d0/
    REAL(KIND=8), DIMENSION(2,temp_mesh%gauss%l_G*temp_mesh%me)                  :: rr_gauss
    INTEGER,      DIMENSION(temp_mesh%gauss%n_w)                               :: j_loc
    !Communicators for Petsc, in space and Fourier------------------------------
#include "petsc/finclude/petsc.h"
    petscerrorcode                   :: ierr
    mpi_comm, DIMENSION(:), POINTER  :: comm_one_d
    mat, DIMENSION(:), POINTER, SAVE :: temp_mat
    vec,                        SAVE :: cb_1, cb_2, cx_1, cx_1_ghost
    ksp, DIMENSION(:), POINTER, SAVE :: temp_ksp
    !------------------------------END OF DECLARATION--------------------------------------

    IF (once) THEN

       once = .false.

       CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)

       !-----CREATE PETSC VECTORS AND GHOSTS-----------------------------------------
       CALL create_my_ghost(temp_mesh,temp_1_la,temp_1_ifrom)
       n = temp_mesh%dom_np
       CALL veccreateghost(comm_one_d(1), n, & 
            petsc_determine, SIZE(temp_1_ifrom), temp_1_ifrom, cx_1, ierr)
       CALL vecghostgetlocalform(cx_1, cx_1_ghost, ierr)
       CALL vecduplicate(cx_1, cb_1, ierr)
       CALL vecduplicate(cx_1, cb_2, ierr)
       !------------------------------------------------------------------------------

       !-------------DIMENSIONS-------------------------------------------------------
       m_max_c = SIZE(list_mode)
       !------------------------------------------------------------------------------

       !---------PREPARE pp_js_D ARRAY FOR TEMPERATURE--------------------------------------
       CALL dirichlet_nodes_parallel(temp_mesh, temp_list_dirichlet_sides, temp_js_d)
       !------------------------------------------------------------------------------

       !--------------------------------------------------------------------------------
       hmoy = 0
       DO m = 1, temp_mesh%dom_me
          hmoy = hmoy + sqrt(sum(temp_mesh%gauss%rj(:,m)))/2
       END DO
       hmoy =  hmoy/temp_mesh%dom_me
       mass0 = 0.d0
       DO i = 1, m_max_c
          mode = list_mode(i)
          IF (mode == 0)  THEN  
             CALL mass_tot(comm_one_d(1),temp_mesh, tempn(:,1,i), mass0)
          ENDIF
       END DO
       !--------------------------------------------------------------------------------

       !-------------ASSEMBLE TEMPERATURE MATRICES--------------------------------------------
       ALLOCATE(temp_mat(m_max_c),temp_ksp(m_max_c))

       DO i = 1, m_max_c
          mode = list_mode(i)

          !---TEMPERATURE MATRIX
          CALL create_local_petsc_matrix(comm_one_d(1), temp_1_la, temp_mat(i), clean=.false.)
          CALL qs_diff_mass_scal_m_variant(temp_mesh, temp_1_la, vol_heat_capacity, temp_diffusivity, &
               1.5/dt, zero, mode, temp_mat(i))

          IF (temp_per%n_bord/=0) THEN
             CALL periodic_matrix_petsc(temp_per%n_bord, temp_per%list, temp_per%perlist, temp_mat(i), temp_1_la)
          END IF

          CALL dirichlet_m_parallel(temp_mat(i),temp_1_la%loc_to_glob(1,temp_js_d))

          CALL init_solver(my_par_cc,temp_ksp(i),temp_mat(i),comm_one_d(1),&
               solver=my_par_cc%solver,precond=my_par_cc%precond)
       ENDDO

    ENDIF
    tps_tot = user_time()
    tps_cumul = 0

    !===Compute rhs by FFT at Gauss points
    tps = user_time()
    CALL smb_ugradc_gauss_fft_par(comm_one_d(2),temp_mesh,list_mode,vol_heat_capacity,chmp_vit,2*tempn-tempn_m1,ff_conv)
    tps = user_time() - tps; tps_cumul=tps_cumul+tps
    !WRITE(*,*) ' Tps fft vitesse', tps
    !------------CONSTRUCTION OF rr_gauss------------------
    index = 0
    DO m = 1, temp_mesh%me
       j_loc = temp_mesh%jj(:,m)
       DO l = 1, temp_mesh%gauss%l_G
          index = index + 1
          rr_gauss(1,index) = sum(temp_mesh%rr(1,j_loc)*temp_mesh%gauss%ww(:,l))
          rr_gauss(2,index) = sum(temp_mesh%rr(2,j_loc)*temp_mesh%gauss%ww(:,l))
       END DO
    END DO
    !------------DEBUT BOUCLE SUR LES MODES----------------
    DO i = 1, m_max_c
       mode = list_mode(i)

       !===RHS temperature
       ff = (2/dt)*tempn(:,1,i) - (1/(2*dt))*tempn_m1(:,1,i)
       CALL qs_00_gauss(temp_mesh, temp_1_la, vol_heat_capacity, ff, &
            -ff_conv(:,1,i) + source_in_temperature(1, rr_gauss, mode, time), cb_1)

       ff = (2/dt)*tempn(:,2,i) - (1/(2*dt))*tempn_m1(:,2,i)
       CALL qs_00_gauss(temp_mesh, temp_1_la, vol_heat_capacity, ff, &
            -ff_conv(:,2,i) + source_in_temperature(2, rr_gauss, mode, time), cb_2)

       !===RHS periodicity
       IF (temp_per%n_bord/=0) THEN
          CALL periodic_rhs_petsc(temp_per%n_bord, temp_per%list, temp_per%perlist, cb_1, temp_1_la)
          CALL periodic_rhs_petsc(temp_per%n_bord, temp_per%list, temp_per%perlist, cb_2, temp_1_la)
       END IF
       !----------------------------------------------------------

       CALL dirichlet_rhs(temp_1_la%loc_to_glob(1,temp_js_d)-1,&
            temperature_exact(1,temp_mesh%rr(:,temp_js_d), mode, time),cb_1)
       CALL dirichlet_rhs(temp_1_la%loc_to_glob(1,temp_js_d)-1,&
            temperature_exact(2,temp_mesh%rr(:,temp_js_d), mode, time),cb_2)

       tps = user_time() - tps; tps_cumul=tps_cumul+tps
       !WRITE(*,*) ' Tps second membre vitesse', tps
       !------------------------------------------------------------------------------------- 

       !--------------------INVERSION DES OPERATEURS--------------
       tps = user_time()
       !Solve system temp_c
       CALL solver(temp_ksp(i),cb_1,cx_1,reinit=.false.,verbose=my_par_cc%verbose)
       CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr) 
       CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
       CALL extract(cx_1_ghost,1,1,temp_1_la,tempn_p1(:,1))

       !Solve system temp_s
       CALL solver(temp_ksp(i),cb_2,cx_1,reinit=.false.,verbose=my_par_cc%verbose)
       CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr) 
       CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
       CALL extract(cx_1_ghost,1,1,temp_1_la,tempn_p1(:,2))
       tps = user_time() - tps; tps_cumul=tps_cumul+tps
       !WRITE(*,*) ' Tps solution des pb de vitesse', tps, 'for mode ', mode
       !------------------------------------------------------------------------------------- 

       !---------------UPDATES-----------------------
       tps = user_time()

       !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010
       IF (mode==0) THEN
          tempn_p1(:,2) = 0.d0
       END IF
       !JLG AR, Dec 18 2008/JLG Bug corrige Jan 23 2010


       tempn_m1(:,:,i) = tempn(:,:,i) 
       tempn(:,:,i) = tempn_p1

       ! Reconstruct dtempndt on Gauss points
       ! WARNING FL (1/2/13) If reactivated, declare l and index
       !tempn_p1 = (3*tempn_p1 - 4*tempn(:,:,i) + tempn_m1(:,:,i))/(2*dt)
       !index = 0
       !DO m = 1, temp_mesh%me
       !   DO l = 1, temp_mesh%gauss%l_G
       !      index  = index + 1
       !      dtempndt(index,1,i) = SUM(tempn_p1(temp_mesh%jj(:,m),1)*temp_mesh%gauss%ww(:,l)) + ff_conv(index,1,i)
       !      dtempndt(index,2,i) = SUM(tempn_p1(temp_mesh%jj(:,m),2)*temp_mesh%gauss%ww(:,l)) + ff_conv(index,2,i)
       !   END DO
       !END DO
       ! WARNING FL (1/2/13) If reactivated, declare l and index

       tps = user_time() - tps; tps_cumul=tps_cumul+tps
       !WRITE(*,*) ' Tps  des updates', tps
       !-------------------------------------------------------------------------------------
    ENDDO

    tps_tot = user_time() - tps_tot
    !WRITE(*,'(A,2(f13.3,2x))') '  Tps boucle en temps Navier_stokes', tps_tot, tps_cumul
    !WRITE(*,*) ' TIME = ', time, '========================================'

  END SUBROUTINE three_level_temperature
  !============================================

  SUBROUTINE smb_ugradc_gauss_fft_par(communicator,mesh,list_mode,heat_capa_in,V_in,c_in,c_out)
    !=================================
    USE gauss_points     
    USE sft_parallele
    USE chaine_caractere
    USE boundary
    IMPLICIT NONE
    TYPE(mesh_type), TARGET                     :: mesh
    INTEGER,      DIMENSION(:),     INTENT(IN)  :: list_mode
    REAL(KIND=8), DIMENSION(:),     INTENT(IN)  :: heat_capa_in
    REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN)  :: v_in, c_in
    REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT) :: c_out
    REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: gradc, w
    REAL(KIND=8), DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: div, cgauss, cint
    INTEGER,      DIMENSION(mesh%gauss%n_w)                  :: j_loc
    REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w)   :: dw_loc     
    INTEGER                                                  ::  m, l , i, mode, index, k
    REAL(KIND=8), DIMENSION(mesh%gauss%n_w,6)   :: vs
    REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2)   :: cs
    INTEGER,      DIMENSION(:,:), POINTER       :: jj
    INTEGER,                      POINTER       :: me
    REAL(KIND=8)   :: ray, tps
    REAL(KIND=8), DIMENSION(3)                  :: temps
    INTEGER                                     :: code, m_max_pad, bloc_size, nb_procs
#include "petsc/finclude/petsc.h"
    mpi_comm       :: communicator

    CALL gauss(mesh)
    jj => mesh%jj
    me => mesh%me

    tps = user_time()
    DO i = 1, SIZE(list_mode)
       mode = list_mode(i)
       index = 0
       DO m = 1, mesh%dom_me
          j_loc = jj(:,m)
          DO k = 1, 6
             vs(:,k) = v_in(j_loc,k,i)
          END DO
          DO k = 1, 2
             cs(:,k) = c_in(j_loc,k,i)
          END DO
          DO l = 1, l_g
             index = index + 1
             dw_loc = dw(:,:,l,m)

             !===Compute radius of Gauss point
             ray = sum(mesh%rr(1,j_loc)*ww(:,l))

             !-----------------vitesse sur les points de Gauss---------------------------
             w(index,1,i) = sum(vs(:,1)*ww(:,l))
             w(index,3,i) = sum(vs(:,3)*ww(:,l))
             w(index,5,i) = sum(vs(:,5)*ww(:,l))

             w(index,2,i) = sum(vs(:,2)*ww(:,l))
             w(index,4,i) = sum(vs(:,4)*ww(:,l))
             w(index,6,i) = sum(vs(:,6)*ww(:,l))

             div(index,1,i) = sum(vs(:,1)*dw_loc(1,:)) + sum(vs(:,1)*ww(:,l))/ray &
                  + (mode/ray)*sum(vs(:,4)*ww(:,l)) +  sum(vs(:,5)*dw_loc(2,:))
             div(index,2,i) = sum(vs(:,2)*dw_loc(1,:)) + sum(vs(:,2)*ww(:,l))/ray &
                  - (mode/ray)*sum(vs(:,3)*ww(:,l)) +  sum(vs(:,6)*dw_loc(2,:))

             !-----------------gradient de c sur les points de Gauss---------------------------
             !coeff sur les cosinus et sinus
             gradc(index,1,i) = sum(cs(:,1)*dw_loc(1,:))
             gradc(index,2,i) = sum(cs(:,2)*dw_loc(1,:))
             gradc(index,3,i) =  mode/ray*sum(cs(:,2)*ww(:,l))
             gradc(index,4,i) = -mode/ray*sum(cs(:,1)*ww(:,l))
             gradc(index,5,i) = sum(cs(:,1)*dw_loc(2,:))
             gradc(index,6,i) = sum(cs(:,2)*dw_loc(2,:))

             gradc(index,:,i) = heat_capa_in(m) * gradc(index,:,i)
             
             cgauss(index,1,i) = sum(cs(:,1)*ww(:,l))
             cgauss(index,2,i) = sum(cs(:,2)*ww(:,l))

             cgauss(index,:,i) = heat_capa_in(m) * cgauss(index,:,i)
             
          ENDDO
       ENDDO
    END DO

    !tps = user_time() - tps
    !WRITE(*,*) ' Tps dans la grande boucle', tps
    !tps = user_time()
    temps = 0

    CALL mpi_comm_size(communicator, nb_procs, code)
    bloc_size = SIZE(gradc,1)/nb_procs+1
    m_max_pad = 3*SIZE(list_mode)*nb_procs/2
    CALL fft_par_dot_prod_dcl(communicator, gradc, w, c_out, nb_procs, bloc_size, m_max_pad, temps)
    bloc_size = SIZE(div,1)/nb_procs+1
    CALL fft_par_prod_dcl(communicator, div, cgauss, cint, nb_procs, bloc_size, m_max_pad, temps)
    c_out = c_out + cint
    tps = user_time() - tps
    !WRITE(*,*) ' Tps dans FFT_PAR_PROD_VECT', tps
    !write(*,*) ' Temps de Comm   ', temps(1)
    !write(*,*) ' Temps de Calc   ', temps(2)
    !write(*,*) ' Temps de Chan   ', temps(3)

  END SUBROUTINE smb_ugradc_gauss_fft_par

  SUBROUTINE mass_tot(communicator,mesh,tempn,RESLT)
    !===========================
    !moyenne
    USE def_type_mesh
    IMPLICIT NONE
    TYPE(mesh_type)                             :: mesh
    REAL(KIND=8), DIMENSION(:)  ,   INTENT(IN)  :: tempn
    REAL(KIND=8)                ,   INTENT(OUT) :: reslt
    REAL(KIND=8)                                :: r_loc, r_out   
    INTEGER ::  m, l , i , ni, code
    INTEGER,      DIMENSION(mesh%gauss%n_w)     :: j_loc
    REAL(KIND=8)   :: ray
#include "petsc/finclude/petsc.h"
    mpi_comm                                    :: communicator
    r_loc = 0.d0

    DO m = 1, mesh%dom_me
       j_loc = mesh%jj(:,m)
       DO l = 1, mesh%gauss%l_G
          ray = 0
          DO ni = 1, mesh%gauss%n_w;  i = j_loc(ni)
             ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
          END DO

          r_loc = r_loc +  sum(tempn(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)

       ENDDO
    ENDDO

    CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
    reslt = r_out 

  END SUBROUTINE mass_tot

END MODULE subroutine_temperature