multi_scale_solve_fe_film Subroutine

  1. mod_ms_film.f90
  2. ms_film
  3. multi_scale_solve_fe_film

public subroutine multi_scale_solve_fe_film(ms_fe_f, ms_mat, bc)

Arguments

Type IntentOptional AttributesName
type(MS_FE_FILM), intent(inout) :: ms_fe_f

multi-scale FE film
type(MS_MAT_SOLV), intent(inout) :: ms_mat

multi-scale solver matrices
real(kind=R8), intent(in), dimension(MAX_NNC):: bc

boundary conditions at the corners

Calls

proc~~multi_scale_solve_fe_film~~CallsGraph proc~multi_scale_solve_fe_film multi_scale_solve_fe_film proc~apply_bc_fe_film_simple apply_bc_FE_film_simple proc~multi_scale_solve_fe_film->proc~apply_bc_fe_film_simple proc~convert_matrice_format convert_matrice_format proc~multi_scale_solve_fe_film->proc~convert_matrice_format omp_get_thread_num omp_get_thread_num proc~multi_scale_solve_fe_film->omp_get_thread_num proc~multi_scale_assembly_fe_film_reynolds multi_scale_assembly_fe_film_reynolds proc~multi_scale_solve_fe_film->proc~multi_scale_assembly_fe_film_reynolds proc~solve_syst solve_syst proc~multi_scale_solve_fe_film->proc~solve_syst proc~solve_fe_film solve_FE_film proc~multi_scale_solve_fe_film->proc~solve_fe_film proc~from_elemental_to_assembled from_elemental_to_assembled proc~convert_matrice_format->proc~from_elemental_to_assembled proc~multi_scale_assembly_fe_film_reynolds->omp_get_thread_num proc~analyse_solver analyse_solver proc~solve_syst->proc~analyse_solver proc~factorize_solver factorize_solver proc~solve_syst->proc~factorize_solver proc~close_solver close_solver proc~solve_syst->proc~close_solver proc~solution_solver solution_solver proc~solve_syst->proc~solution_solver proc~freefact_solver freefact_solver proc~solve_syst->proc~freefact_solver proc~init_solver init_solver proc~solve_syst->proc~init_solver proc~solve_fe_film->proc~apply_bc_fe_film_simple proc~solve_fe_film->proc~convert_matrice_format proc~solve_fe_film->omp_get_thread_num proc~solve_fe_film->proc~solve_syst proc~apply_bc_fe_film apply_bc_FE_film proc~solve_fe_film->proc~apply_bc_fe_film proc~assembly_fe_film_reynolds assembly_FE_film_reynolds proc~solve_fe_film->proc~assembly_fe_film_reynolds ma48_analyse ma48_analyse proc~analyse_solver->ma48_analyse proc~prep_superlu prep_superlu proc~analyse_solver->proc~prep_superlu dmumps dmumps proc~analyse_solver->dmumps proc~s_umfpack_di_symbolic s_umfpack_di_symbolic proc~analyse_solver->proc~s_umfpack_di_symbolic proc~s_umfpack_di_numeric s_umfpack_di_numeric proc~factorize_solver->proc~s_umfpack_di_numeric proc~umfpack_di_free_numeric umfpack_di_free_numeric proc~factorize_solver->proc~umfpack_di_free_numeric proc~factorize_solver->dmumps ma48_factorize ma48_factorize proc~factorize_solver->ma48_factorize proc~umfpack_di_free_symbolic umfpack_di_free_symbolic proc~close_solver->proc~umfpack_di_free_symbolic proc~close_solver->proc~umfpack_di_free_numeric proc~close_solver->dmumps mpi_finalize mpi_finalize proc~close_solver->mpi_finalize proc~close_superlu close_superlu proc~close_solver->proc~close_superlu proc~umfpack_di_report_info umfpack_di_report_info proc~close_solver->proc~umfpack_di_report_info ma48_finalize ma48_finalize proc~close_solver->ma48_finalize proc~solv_superlu solv_superlu proc~solution_solver->proc~solv_superlu ma48_solve ma48_solve proc~solution_solver->ma48_solve proc~solution_solver->dmumps proc~solution_solver->mpi_finalize proc~s_umfpack_di_solve s_umfpack_di_solve proc~solution_solver->proc~s_umfpack_di_solve proc~freefact_solver->proc~umfpack_di_free_numeric proc~free_superlu free_superlu proc~freefact_solver->proc~free_superlu proc~elementary_assembly_fe_film_reynolds elementary_assembly_FE_film_reynolds proc~assembly_fe_film_reynolds->proc~elementary_assembly_fe_film_reynolds proc~assemble_in_mat_sol assemble_in_mat_sol proc~assembly_fe_film_reynolds->proc~assemble_in_mat_sol mpi_init mpi_init proc~init_solver->mpi_init ma48_initialize ma48_initialize proc~init_solver->ma48_initialize proc~init_solver->dmumps proc~umfpack_di_report_control umfpack_di_report_control proc~init_solver->proc~umfpack_di_report_control proc~init_superlu init_superlu proc~init_solver->proc~init_superlu proc~init_solver->mpi_finalize proc~umfpack_di_defaults umfpack_di_defaults proc~init_solver->proc~umfpack_di_defaults proc~umfpack_di_numeric umfpack_di_numeric proc~s_umfpack_di_numeric->proc~umfpack_di_numeric proc~solv_superlu->proc~prep_superlu interface~statfree StatFree proc~solv_superlu->interface~statfree interface~destroy_supernode_matrix Destroy_SuperNode_Matrix proc~solv_superlu->interface~destroy_supernode_matrix interface~destroy_compcol_matrix Destroy_CompCol_Matrix proc~solv_superlu->interface~destroy_compcol_matrix interface~dgssvx dgssvx proc~solv_superlu->interface~dgssvx interface~statinit StatInit proc~solv_superlu->interface~statinit interface~dcreate_dense_matrix dCreate_Dense_Matrix proc~prep_superlu->interface~dcreate_dense_matrix interface~dcreate_compcol_matrix dCreate_CompCol_Matrix proc~prep_superlu->interface~dcreate_compcol_matrix interface~c_umfpack_di_free_symbolic c_umfpack_di_free_symbolic proc~umfpack_di_free_symbolic->interface~c_umfpack_di_free_symbolic interface~c_umfpack_di_free_numeric c_umfpack_di_free_numeric proc~umfpack_di_free_numeric->interface~c_umfpack_di_free_numeric interface~c_umfpack_di_report_control c_umfpack_di_report_control proc~umfpack_di_report_control->interface~c_umfpack_di_report_control interface~set_default_options set_default_options proc~init_superlu->interface~set_default_options proc~compute_prc_tables_reynolds_supg compute_prc_tables_reynolds_supg proc~elementary_assembly_fe_film_reynolds->proc~compute_prc_tables_reynolds_supg proc~umfpack_di_solve umfpack_di_solve proc~s_umfpack_di_solve->proc~umfpack_di_solve interface~destroy_dense_matrix Destroy_Dense_Matrix proc~close_superlu->interface~destroy_dense_matrix proc~close_superlu->interface~destroy_compcol_matrix proc~umfpack_di_symbolic umfpack_di_symbolic proc~s_umfpack_di_symbolic->proc~umfpack_di_symbolic interface~c_umfpack_di_report_info c_umfpack_di_report_info proc~umfpack_di_report_info->interface~c_umfpack_di_report_info interface~c_umfpack_di_defaults c_umfpack_di_defaults proc~umfpack_di_defaults->interface~c_umfpack_di_defaults proc~ni4_up_2d ni4_up_2d proc~compute_prc_tables_reynolds_supg->proc~ni4_up_2d proc~length_width_elem length_width_elem proc~compute_prc_tables_reynolds_supg->proc~length_width_elem proc~dj4 dj4 proc~compute_prc_tables_reynolds_supg->proc~dj4 interface~c_umfpack_di_symbolic c_umfpack_di_symbolic proc~umfpack_di_symbolic->interface~c_umfpack_di_symbolic interface~c_umfpack_di_numeric c_umfpack_di_numeric proc~umfpack_di_numeric->interface~c_umfpack_di_numeric interface~c_umfpack_di_solve c_umfpack_di_solve proc~umfpack_di_solve->interface~c_umfpack_di_solve proc~ni4_up_1d ni4_up_1d proc~ni4_up_2d->proc~ni4_up_1d

Called by

proc~~multi_scale_solve_fe_film~~CalledByGraph proc~multi_scale_solve_fe_film multi_scale_solve_fe_film proc~solve_ms_prob solve_ms_prob proc~solve_ms_prob->proc~multi_scale_solve_fe_film proc~test_slider_ms test_slider_ms proc~test_slider_ms->proc~solve_ms_prob proc~test_rough_ms test_rough_ms proc~test_rough_ms->proc~solve_ms_prob proc~run_test run_test proc~run_test->proc~test_slider_ms proc~run_test->proc~test_rough_ms program~main main program~main->proc~run_test

Contents

Source Code

multi_scale_solve_fe_film

Source Code

   subroutine multi_scale_solve_fe_film(ms_fe_f, ms_mat, bc)
   implicit none
   type(MS_FE_FILM),  intent(inout)                     :: ms_fe_f   !! *multi-scale FE film*
   type(MS_MAT_SOLV), intent(inout)                     :: ms_mat    !! *multi-scale solver matrices*
   real(kind=R8),     intent(in   ), dimension(MAX_NNC) :: bc        !! *boundary conditions at the corners*

      logical(kind=I4) :: decomp
      integer(kind=I4) :: i, e, it, fluid
      logical(kind=I4) :: conv, ts_init
      real(kind=R8)    :: error, relax
      integer(kind=I4) :: ass_c
      integer(kind=I4), dimension(2) :: compt

      relax = ms_fe_f%ts_fe_f%num_p%relax
      fluid = ms_fe_f%ts_fe_f%data_f%fl%fluid_type

      VERBOSE = nint(VERBOSE/10.)

      ! open mp instructions (parallele computation)
      if (VERBOSE >= 1) write(OPU,*) 'nb_threads_max', NB_THREADS_MAX
      if (NB_THREADS_MAX <= 0) NB_THREADS_MAX = 1

      !$ call omp_set_num_threads(NB_THREADS_MAX)
      if (VERBOSE >= 1) write(OPU,*) 'nb_threads_used', NB_THREADS_MAX


      ! check of matrices allocation
      allocate(ms_mat%bs_mat(ms_fe_f%ts_fe_f%m%ne))
      ms_mat%bs_mat(:)%slv_t = SOLVER_BS
      ms_mat%bs_mat(:)%first = .true.

      ! solution parameters
      ass_c  = ASS
      decomp = (ass_c == ASS)

      ! update fluid properties
      do i = 1, ms_fe_f%ts_fe_f%m%n
         ms_fe_f%ts_fe_f%vn(i,RHO_N) = ms_fe_f%ts_fe_f%data_f%fl%rho( ms_fe_f%ts_fe_f%vn(i, P_N), &
                                                                      ms_fe_f%ts_fe_f%vn(i, T_N) )

         ms_fe_f%ts_fe_f%vn(i,DRHODP_N) = ms_fe_f%ts_fe_f%data_f%fl%drhodp( ms_fe_f%ts_fe_f%vn(i, P_N), &
                                                                            ms_fe_f%ts_fe_f%vn(i, T_N) )
      enddo
      if (VERBOSE >= 2) write(OPU,*) 'fluid properties updated'

      ms_mat%ts_mat%slv_t = SOLVER_TS
      ms_mat%ts_mat%first = .true.
      ms_mat%ts_mat%nn = ms_fe_f%ts_fe_f%m%n
      ms_mat%ts_mat%ne = ms_fe_f%ts_fe_f%m%ne

      ! matrices allocation
      compt(:) = 0
      do e = 1, ms_fe_f%ts_fe_f%m%ne
         compt(1) = compt(1) +ms_fe_f%ts_fe_f%m%el_n(e)       ! for each element, the number of lines is added
         do i = 1, ms_fe_f%ts_fe_f%m%el_t(e)
            compt(2) = compt(2) + ms_fe_f%ts_fe_f%m%el_t(e)   ! " ", for each node, the number of contributions is added
         enddo
      enddo
      ms_mat%ts_mat%nvar = compt(1)
      ms_mat%ts_mat%nt   = compt(2)
      allocate( ms_mat%ts_mat%eltvar(ms_mat%ts_mat%nvar  ) )
      allocate( ms_mat%ts_mat%a_elt( ms_mat%ts_mat%nt    ) )
      allocate( ms_mat%ts_mat%eltptr(ms_mat%ts_mat%ne +1 ) )

      ! check of precomputed tables allocation
      if (.not.allocated(ms_fe_f%ts_fe_f%prc%vcal)) call init_prc_tab(ms_fe_f%ts_fe_f)

      ! convergence is false at zero iteration
      conv = .false.
      it   = 0

      ts_init = .false.
      if (ts_init) then
         call solve_FE_film(fe_f = ms_fe_f%ts_fe_f,   &
                             mat = ms_mat%ts_mat,     &
                              bc = bc,                &
                        flag_ass = .true.)
      else
         ! apply boundary conditions
         call apply_bc_FE_film_simple(fe_f = ms_fe_f%ts_fe_f,   &
                                        bc = bc)
         if (VERBOSE >= 2) write(OPU,*) 'bc applied'
         call solve_syst(mat = ms_mat%ts_mat, &
                        step = 'ini')
         if (VERBOSE >= 2) write(OPU,*) 'TS solver initialized'
      endif

      if ( sum(ms_fe_f%bs_fe_f(1)%m%ed(:)%n) > MAX_NBS ) stop 'MAX_NBS under estimated'

      ! solution loop
      do
         if (VERBOSE >= 3) write(OPU,*) "loop MS *******************", it

         if (conv) then
            do e = 1, ms_fe_f%ts_fe_f%m%ne
               call solve_syst(mat = ms_mat%bs_mat(e), & !
                              step = 'end')
               if (VERBOSE >= 3) write(OPU,*) '   Matrix BS released, thread ', omp_get_thread_num()
            enddo
            exit
         endif

         if (SMOOTH_MS) call smooth_ms_fe_f(ms_fe_f, code=P_N, nodal=.true.)

         if (BC_SPLINE) call interp_ts_bs_splin(ms_fe_f)

         ! assembly of the system
         call multi_scale_assembly_fe_film_reynolds(ms_fe_f = ms_fe_f, &
                                                    ms_mat  = ms_mat,  &
                                                      ass_c = ass_c)
         if (VERBOSE >= 3) write(OPU,*) 'ms reynolds assembled'
         if (VERBOSE >= 3) write(OPU,*) 'ass_c', ass_c, 'first', ms_mat%ts_mat%first

!~          if (ass_c == ASS) then
            ! some stuff can be saved here, provided the reloading of jptr, irow, ... (instead of convert_matrice_format)
            call convert_matrice_format(mat = ms_mat%ts_mat)
            if (VERBOSE >= 3) write(OPU,*) 'Matrix TS formated, thread ', omp_get_thread_num()
!~          endif

         if (ms_mat%ts_mat%first) then
            call solve_syst(mat = ms_mat%ts_mat, &
                           step = 'ana')
            ms_mat%ts_mat%first = .false.
            if (VERBOSE >= 3) write(OPU,*) 'Matrix TS analyzed, thread ', omp_get_thread_num()
         endif

         ! solution of the system
!~          if (ass_c == ASS) then
            call solve_syst(mat = ms_mat%ts_mat, &
                           step = 'fac')
            if (VERBOSE >= 3) write(OPU,*) 'Matrix TS factorized, thread ', omp_get_thread_num()
!~          endif
         call solve_syst(mat = ms_mat%ts_mat, &
                        step = 'sol')
         if (VERBOSE >= 3) write(OPU,*) 'System TS solved, thread ', omp_get_thread_num()

         !if (ass_c == ASS) then
            call solve_syst(mat = ms_mat%ts_mat, &
                           step = 'fre')
            if (VERBOSE >= 3) write(OPU,*) 'Matrix factors freed, thread ', omp_get_thread_num()
         !endif

         ! error computation
         error = (sum(ms_mat%ts_mat%x ** 2) / sum(ms_fe_f%ts_fe_f%vn(:, P_N) ** 2)) ** (0.5)
         it = it + 1
         if (VERBOSE >= 1) write(OPU,*) 'Iteration ', it, 'error TS', error

         ! convergence check
         if (error <= ms_fe_f%ts_fe_f%num_p%eps) conv = .true.

         ! update of variables
         if (fluid == MIXT) then
            do i = 1, ms_fe_f%ts_fe_f%m%n
               if (ms_mat%ts_mat%x(i) < 0.) then
                  ms_fe_f%ts_fe_f%vn(i, RHO_N) = ms_fe_f%ts_fe_f%vn(i, RHO_N) + ms_fe_f%ts_fe_f%vn(i, DRHODP_N) * ms_mat%ts_mat%x(i) * relax
                  if (ms_fe_f%ts_fe_f%vn(i, RHO_N) < 0.) then
                      ms_fe_f%ts_fe_f%vn(i,   P_N) = ms_fe_f%ts_fe_f%data_f%fl%p_0/100
                  else
                      ms_fe_f%ts_fe_f%vn(i,   P_N) = ms_fe_f%ts_fe_f%data_f%fl%pres( ms_fe_f%ts_fe_f%vn(i, RHO_N), & !
                                                                                     ms_fe_f%ts_fe_f%vn(i,   T_N) )
                  endif
               else
                  ms_fe_f%ts_fe_f%vn(i, P_N) = ms_fe_f%ts_fe_f%vn(i, P_N) + ms_mat%ts_mat%x(i) * relax
               endif
            enddo
         else
            ms_fe_f%ts_fe_f%vn(:, P_N) = ms_fe_f%ts_fe_f%vn(:, P_N) + ms_mat%ts_mat%x * relax
         endif

         ! check pressure
         if ( fluid == GP ) then
            if (minval(ms_fe_f%ts_fe_f%vn(:, P_N)) < 0._R8) write(OPU,*) 'MS p negative'
            where (ms_fe_f%ts_fe_f%vn(:, P_N) < 0._R8) ms_fe_f%ts_fe_f%vn(:, P_N) = ms_fe_f%ts_fe_f%data_f%fl%p_0 / 1.e2_R8
         endif

         ! update fluid properties
         do i = 1, ms_fe_f%ts_fe_f%m%n
            ms_fe_f%ts_fe_f%vn(i, RHO_N) = ms_fe_f%ts_fe_f%data_f%fl%rho( ms_fe_f%ts_fe_f%vn(i, P_N), &
                                                                          ms_fe_f%ts_fe_f%vn(i, T_N) )

            ms_fe_f%ts_fe_f%vn(i, DRHODP_N) = ms_fe_f%ts_fe_f%data_f%fl%drhodp( ms_fe_f%ts_fe_f%vn(i, P_N), &
                                                                                ms_fe_f%ts_fe_f%vn(i, T_N) )
         enddo

         if (it >= ms_fe_f%ts_fe_f%num_p%it_max) then
            conv = .true.
            write(OPU,*) 'maximum number of iteration reached before convergence'
         endif


      enddo

      call solve_syst(mat = ms_mat%ts_mat, & !
                     step = 'end')
      if (VERBOSE >= 2) write(OPU,*) 'Matrix TS released, thread ', omp_get_thread_num()

   return
   endsubroutine multi_scale_solve_fe_film