program STOKES_1p3

!     Escoamento Bidimensional Isotérmico e
!     Permanente de Fluido Incompressível

!     Versão atual    1.3 ( 6 Out 05)

!     Versão          1.2 (29 Jul 05)
!     Versão          1.1 (30 Abr 99)
!     Versão original 1.0 (25 Mar 99)

!     última alteração =   11 Out 05

!     Autor: Carlos Henrique Marchi
!            DEMEC/UFPR, Curitiba, PR


! MODELO MATEMÁTICO (equações diferenciais):
! -----------------------------------------

!   massa: du/dx + dv/dy = 0

!   QMLx:  d(u2)/dx + d(uv)/dy = -dp/dx + mi*laplaciano(u)

!   QMLy:  d(uv)/dx + d(v2)/dy = -dp/dy + mi*laplaciano(v) - S

!   x, y = coordenadas espaciais
!   u, v = velocidades nas direções x e y
!   S    = termo fonte; nulo no problema de Ghia
!   Re   = Ut * L / mi (número de Reynolds)
!   Ut   = velocidade máxima da tampa da cavidade (1 m/s)
!   Cavidade quadrada (L=1 por H=1 m)
!   ro   = massa específica (1 kg/m3)
!   Condição de contorno de Dirichlet para u e v
!   Tipos de problemas pré-definidos: 
!     1 = Shih: u na tampa é uma função de x; u e v nulos nas outras paredes
!     2 = Ghia: u na tampa é constante = 1; u e v nulos nas outras paredes


! MODELO NUMÉRICO (resumo):
! ------------------------

!   Incógnitas: u, v, p
!   Condições de contorno de Dirichlet para u e v, e p, extrapolado;
!     aplicadas com volumes fictícios
!   Método numérico: volumes finitos
!   Funções de interpolação para u e v: 
!     CDS para os termos difusivos
!     UDS ou CDS (com correção adiada) para os termos advectivos
!   Malha uniforme em cada direção
!   Método de acoplamento pressão-velocidade: SIMPLEC
!   Solução segregada (seqüencial) das equações diferenciais
!   Solvers: Jacobi, Gauss-Seidel, SOR, ADI e MSI-5
!   Precisão: dupla
!   Tempo é usado para avanço do processo iterativo
!   Arranjo co-localizado de variáveis (média NHT/1994)
!   Linguagem FORTRAN 95
!   Aplicativo usado: Compaq Visual Fortran 6.6
!   Tipo de projeto: QuickWin Application
!   Variáveis de interesse:
!     Fluxo de massa que escoa na cavidade: obtido com a função de corrente
!        baseada na velocidade média u na face leste de cada volume
!     Fluxo de massa que escoa em y=0.5 e x[0;0.5]: obtido com a integração 
!        pela regra do retângulo da velocidade média v na face norte de cada volume
!     Força da placa sobre o fluido: obtida com integração pela regra do retângulo 
!        da derivada du/dy em y=1, obtida com velocidade nodal u e com UDS e UDS-2
!     Derivada du/dy em x=0.5 e y=1: obtida com velocidade nodal u 
!        e com UDS e UDS-2; válido só para número de volumes ímpares em x
!     Derivada du/dy em x=0.5 e y=1: obtida com velocidade média u 
!        e com UDS e UDS-2; válido só para número de volumes pares em x;
!        cálculos desativados
!     Pressão no centro do domínio: obtido por média aritmética de p nos
!        quatro nós vizinhos
!     Perfil da velocidade u em x = 0.5: obtido por média aritmética das
!        velocidades na face leste nos dois volumes vizinhos
!     Perfil da velocidade v em y = 0.5: obtido por média aritmética das
!        velocidades na face norte nos dois volumes vizinhos
!   Estimativa inicial da solução numérica é igual a:
!     Solução analítica do problema de Shih ou Zero
!   Expressão genérica dos sistemas de equações discretizados:
!     a(NP,1)*T(NPS) + a(NP,2)*T(NPW) + a(NP,3)*T(NP) + a(NP,4)*T(NPE) 
!                    + a(NP,5)*T(NPN) = b(NP)
!     onde T representa u, v ou plinha


! ARQUIVOS envolvidos no programa:
! -------------------------------

!   Executável:
!     Stokes_1p3.exe

!   Dados:
!     Stokes_1p3.ent
!     u_Shih.gnu = dados para gerar gráfico dos perfis de u no problema de Shih
!     v_Shih.gnu = dados para gerar gráfico dos perfis de v no problema de Shih
!     u_Ghia.gnu = dados para gerar gráfico dos perfis de u no problema de Ghia
!     v_Ghia.gnu = dados para gerar gráfico dos perfis de v no problema de Ghia
!     R.gnu      = dados para gerar gráfico dos resíduo das equações ao longo das iterações
!     Ghia_Re_100_U.dat   = perfil de u no problema de Ghia para Re = 100
!     Ghia_Re_1000_U.dat  = perfil de u no problema de Ghia para Re = 1000
!     Ghia_Re_10000_U.dat = perfil de u no problema de Ghia para Re = 10000
!     Ghia_Re_100_V.dat   = perfil de v no problema de Ghia para Re = 100
!     Ghia_Re_1000_V.dat  = perfil de v no problema de Ghia para Re = 1000
!     Ghia_Re_10000_V.dat = perfil de v no problema de Ghia para Re = 10000

!   Saídas:
!     "caso"     = saída principal do programa
!     "caso".itn = variáveis de interesse ao longo das iterações
!     r.dat      = resíduo das equações ao longo das iterações
!     ucex.dat   = perfil de u em x=1/2 (exato)
!     vcex.dat   = perfil de v em y=1/2 (exato).
!     u.dat      = perfil de u em x=1/2 (numérico)
!     v.dat      = perfil de v em y=1/2 (numérico)

!   Fontes:
!     Stokes_1p3.f90 = programa principal
!     coef.f90       = calcula coeficientes dos sistemas de equações
!     dados.f90      = lê e lista os dados do programa
!     result.f90     = resolve equações e gera listagem dos resultados
!     solvers.f90    = solvers disponíveis
!     varia.f90      = define todas as variáveis globais do programa
!
!   Outros arquivos:
!     Notepad.exe    = editor de arquivos
!     Wgnuplot.exe   = programa gerador de gráfico (versão 3.5)

! -----------------------------------------------

use resultados

! -----------------------------------------------

implicit none

!-------------------------------------------------

interface
  logical(4) function initialsettings
  end function
end interface

  name=setwindowtext(gethwndqq(qwin$framewindow),  &
       "STOKES 1.3 (Escoamento Bidimensional de Fluido Incompressível)"C)

! -----------------------------------------------

  call le_dados

  open(10,file=caso)

  call data_hora

  call mostra_dados

! -----------------------------------------------

  call solucao_numerica

  close (10)

  note_caso = 'notepad '//caso
  dos = system(note_caso) ! lista arquivo de resultados

! -----------------------------------------------

fim = setexitqq(qwin$exitnopersist) 

end

!--------------------------------------------------------

logical(4) function initialsettings()

  use msflib

  record /qwinfo/ qwi
  qwi.x = 100
  qwi.y =  20
  qwi.w = 600         ! 400  600
  qwi.h = 400         ! 200  400
  qwi.type = qwin$set
  i4 = setwsizeqq(qwin$framewindow,qwi)
  qwi.type = qwin$min
  i4 = setwsizeqq(qwin$framewindow,qwi)

  initialsettings = .true.

end function initialsettings

!--------------------------------------------------------
