compute_rot_h.f90

Go to the documentation of this file.

  SUBROUTINE compute_rot_h(field_in, field_out)
    USE def_type_mesh
    USE my_util
    USE fem_m_axi
    USE st_matrix
    USE solve_petsc
    IMPLICIT NONE
    REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN)    :: field_in
    REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: field_out

    INTEGER, POINTER, DIMENSION(:)                :: ifrom
    TYPE(solver_param)                            :: param_rot_h
    REAL(KIND=8), DIMENSION(SIZE(field_out,1), SIZE(field_out,2)) :: smb_rot_h
    REAL(KIND=8), DIMENSION(6)                  :: rot_h_loc
    INTEGER, DIMENSION(mesh_rot_h%gauss%n_w) :: i_loc
    INTEGER               :: i, mode, m, l, j, j_loc, k, nj
    REAL(KIND=8)          :: ray, drdthetadz
#include "petsc/finclude/petsc.h"
    petscerrorcode                              :: ierr


    param_rot_h%it_max=1
    param_rot_h%rel_tol=1.d-6
    param_rot_h%abs_tol=1.d-10
    param_rot_h%solver='GMRES'
    param_rot_h%precond='MUMPS'
    param_rot_h%verbose=.false.

    IF (if_first_rot_h) THEN
       if_first_rot_h=.false.
       
       CALL create_my_ghost(mesh_rot_h, vizu_rot_h_la,ifrom)
       CALL veccreateghost(vizu_grad_curl_comm(1), mesh_rot_h%dom_np, & 
            petsc_determine, SIZE(ifrom), ifrom, vx_phi, ierr)
       CALL vecghostgetlocalform(vx_phi, vx_phi_ghost, ierr)
       CALL vecduplicate(vx_phi, vb_phi, ierr)

       CALL create_local_petsc_matrix(vizu_grad_curl_comm(1),vizu_rot_h_la , mat_rot_h, clean=.true.)
       CALL qs_diff_mass_scal_m (mesh_rot_h, vizu_rot_h_la, 0.d0, 1.d0, 0.d0, 0, mat_rot_h)
       CALL init_solver(param_rot_h,ksp_rot_h,mat_rot_h,vizu_grad_curl_comm(1),&
            solver='GMRES',precond='MUMPS')
       CALL matdestroy(mat_rot_h,ierr)

    END IF
    
    smb_rot_h = 0
    field_out = 0.d0
    DO i = 1, SIZE(vizu_list_mode)
       mode = vizu_list_mode(i)
       smb_rot_h=0.d0
       DO m = 1, mesh_rot_h%me
          DO l = 1, mesh_rot_h%gauss%l_G
             ray = sum(mesh_rot_h%rr(1,mesh_rot_h%jj(:,m))*mesh_rot_h%gauss%ww(:,l))
             drdthetadz = mesh_rot_h%gauss%rj(l,m)
             rot_h_loc = 0.d0
             DO nj = 1,mesh_rot_h%gauss%n_w
                j_loc = mesh_rot_h%jj(nj,m)
! CN attention, here is ray*rot_h_loc
!             RotH(index,1,i) = mode/ray*H_gauss(index,6,i) &
!                  -SUM(Hs(:,3)*dw_loc(2,:))
!             RotH(index,2,i) =-mode/ray*H_gauss(index,5,i) &
!                  -SUM(Hs(:,4)*dw_loc(2,:))
!             RotH(index,3,i) =         SUM(Hs(:,1)*dw_loc(2,:)) &
!                  -SUM(Hs(:,5)*dw_loc(1,:))
!             RotH(index,4,i) =          SUM(Hs(:,2)*dw_loc(2,:)) &
!                  -SUM(Hs(:,6)*dw_loc(1,:))
!             RotH(index,5,i) =    1/ray*H_gauss(index,3,i) &
!                  +SUM(Hs(:,3)*dw_loc(1,:)) &
!                  -mode/ray*H_gauss(index,2,i)
!             RotH(index,6,i) =    1/ray*H_gauss(index,4,i) &
!                  +SUM(Hs(:,4)*dw_loc(1,:))&
!                  +mode/ray*H_gauss(index,1,i)
!
                rot_h_loc(1) = rot_h_loc(1) + mode*field_in(j_loc,6,i) &
                              -field_in(j_loc,3,i)*mesh_rot_h%gauss%dw(2,nj,l,m)*ray
                rot_h_loc(2) = rot_h_loc(2) - mode*field_in(j_loc,5,i) &
                              -field_in(j_loc,4,i)*mesh_rot_h%gauss%dw(2,nj,l,m)*ray
                rot_h_loc(3) = rot_h_loc(3) + field_in(j_loc,1,i)*mesh_rot_h%gauss%dw(2,nj,l,m)*ray &
                              -field_in(j_loc,5,i)*mesh_rot_h%gauss%dw(1,nj,l,m)*ray
                rot_h_loc(4) = rot_h_loc(4) + field_in(j_loc,2,i)*mesh_rot_h%gauss%dw(2,nj,l,m)*ray &
                              -field_in(j_loc,6,i)*mesh_rot_h%gauss%dw(1,nj,l,m)*ray
                rot_h_loc(5) = rot_h_loc(5) + field_in(j_loc,3,i) &
                              +field_in(j_loc,3,i)*mesh_rot_h%gauss%dw(1,nj,l,m)*ray &
                              - mode*field_in(j_loc,2,i)

                rot_h_loc(6) = rot_h_loc(6) + field_in(j_loc,4,i) &
                              +field_in(j_loc,4,i)*mesh_rot_h%gauss%dw(1,nj,l,m)*ray &
                              + mode*field_in(j_loc,1,i)
             END DO
             i_loc = mesh_rot_h%jj(:,m)

             DO k = 1, 6
                smb_rot_h(i_loc,k) = smb_rot_h(i_loc,k) + rot_h_loc(k)* &
                     mesh_rot_h%gauss%ww(:,l)*drdthetadz
             END DO
          END DO
       END DO
 
       DO k = 1, 6
          CALL veczeroentries(vb_phi, ierr)
          CALL vecsetvalues(vb_phi, mesh_rot_h%np, vizu_rot_h_la%loc_to_glob(1,:)-1, smb_rot_h(:,k), add_values, ierr)
          CALL vecassemblybegin(vb_phi,ierr)
          CALL vecassemblyend(vb_phi,ierr)

          CALL solver(ksp_rot_h,vb_phi,vx_phi,reinit=.false.,verbose=.false.)

          CALL vecghostupdatebegin(vx_phi,insert_values,scatter_forward,ierr) 
          CALL vecghostupdateend(vx_phi,insert_values,scatter_forward,ierr)
          IF (mesh_rot_h%me/=0) THEN
             CALL extract(vx_phi_ghost,1,1,vizu_rot_h_la,field_out(:,k,i))
          END IF
       END DO
       
       IF (mode == 0) THEN
          field_out(:,2,i) = 0.d0
          field_out(:,4,i) = 0.d0
          field_out(:,6,i) = 0.d0
       END IF

    END DO
    

  END SUBROUTINE compute_rot_h