grad T forçado ok

anotações.txt

@@ -48,3 +48,21 @@ Ideias:
 
 Variar densidade (número)
 Variar Temperatura
+
+In [1]: m_C = 12
+
+In [2]: m_Ar = 40
+
+In [3]: fe = (.74*(.335+.127)) # ou 0.64*
+
+In [4]: rs = 10 
+
+In [5]: rC = 0.07
+
+In [6]: sigma_Ar = 0.341
+
+In [7]: nC = (sigma_Ar*rs**2*fe/(rC**2) # numero de atomos na região
+
+In [8]: m_Cnd = nC*m_C/m_Ar # massa adimensional do carbono
+
+In [9]: print(nC)

debug_lennard.sh

@@ -1,20 +1,28 @@
 #!/bin/bash
-read -p "[d]ebug with mpiifort or [c]ompile with OpenMPI+gfortran: " option
+read -p "[d]ebug with mpiifort, compile with mpi[i]fort or [c]ompile with OpenMPI+gfortran: " option
 
-if [option = 'c']
+if [ $option == 'c' ] 
 then
     echo "compiling with mpifort"
     /usr/lib64/mpi/gcc/openmpi/bin/mpifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -o lennard.out
     echo "output: lennard.out"
     read -p "Run lennard.out? [y/n]" option2
-    if [option2 = 'y']
+    if [ $option2 = 'y' ]
     then
         /usr/lib64/mpi/gcc/openmpi/bin/mpirun -n 1 ./lennard.out
     fi
 else
-    echo "debuging with intel mpiifort"
-    mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -O0 -shared-intel -nocheck -fno-omit-frame-pointer -fasynchronous-unwind-tables -fexceptions -o lennard.out
-    echo "output: lennard.out"
+    if [ $option == "d" ]
+    then
+        echo "debuging with intel mpiifort"
+        mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -O0 -shared-intel -nocheck -fno-omit-frame-pointer -fasynchronous-unwind-tables -fexceptions -o lennard.out
+        echo "output: lennard.out"
+    else
+        echo "compiling with intel mpiifort without options"
+        mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -o lennard.out
+        echo "output: lennard"
+    fi
+
 fi
 
 

lennard.f90

@@ -466,7 +466,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu)
         ! print*, "FIM"
     end subroutine comp_pot    
 
-    function comp_Kglobal(malha,domx,domy,propriedade,np,id,t) result(K)
+    function comp_Kglobal(malha,domx,domy,propriedade,np,id,t,ignore_mpi) result(K)
         use linkedlist
         use mod1
         use data
@@ -480,8 +480,16 @@ function comp_Kglobal(malha,domx,domy,propriedade,np,id,t) result(K)
         real(dp) :: kb = 1.38064852E-23,K, auxres(8), nump, aux(2)
         type(list_t), pointer :: node
         integer ( kind = 4 ) :: np, ierr, id
+        logical, intent(in), OPTIONAL :: ignore_mpi
+        logical :: i_m
         real(dp), intent(in) :: t
 
+        if (present(ignore_mpi)) then 
+            i_m = ignore_mpi 
+        else
+            i_m = .false.
+        end if
+
         K = 0; nump = 0;
 
         do i = domy(1),domy(2)
@@ -503,7 +511,7 @@ function comp_Kglobal(malha,domx,domy,propriedade,np,id,t) result(K)
         ! então o processo é paralelo, vamos juntar tudo 
         aux = [K, nump]
 
-        if (np > 1) then
+        if (np > 1 .and. (.not. i_m)) then
             call MPI_GATHER(aux, 2, MPI_DOUBLE_PRECISION, auxres, 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
             if (id == 0) then
                 K = auxres(1)+auxres(3)+auxres(5)+auxres(7)
@@ -517,7 +525,7 @@ function comp_Kglobal(malha,domx,domy,propriedade,np,id,t) result(K)
     end function comp_Kglobal
 
     !Corrige energia cinética/configura temepratura global
-    subroutine corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,tt)
+    subroutine corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,tt, ignore_mpi)
         use linkedlist
         use mod1
         use data
@@ -533,9 +541,14 @@ subroutine corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,tt)
         type(list_t), pointer :: node
         integer ( kind = 4 ) :: np, id
         real(dp), intent(in) :: tt
+        logical, INTENT(IN), OPTIONAL :: ignore_mpi
         !calcula temperatura atual
         ! T = (2/(Nf*kb))*Ekin com Nf = número de graus de liberadade 
-        T = comp_Kglobal(malha,domx,domy,propriedade,np,id,tt)
+        if (present(ignore_mpi)) then 
+            T = comp_Kglobal(malha,domx,domy,propriedade,np,id,tt,ignore_mpi)
+        else
+            T = comp_Kglobal(malha,domx,domy,propriedade,np,id,tt)
+        end if
 
         beta = sqrt(Td/T) ! aqui o T é o Beta
         do i = domy(1),domy(2)
@@ -623,6 +636,7 @@ subroutine corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedad
         type(list_t), pointer :: node
         integer ( kind = 4 ) :: np, id
         real(dp), intent(in) :: tt
+
         !calcula temperatura atual
         ! T = (2/(Nf*kb))*Ekin com Nf = número de graus de liberadade 
         T_hot =  comp_K(malha,domx,domy,hot_cells,propriedade,np,id,tt)
@@ -4128,7 +4142,7 @@ program main
     implicit none
 !    Variáveis
     integer :: N,Ntype,i=1, nimpre,j = 1, ii, quant = 0,mesh(2), cont = 1, aux1 = 0,cont2 = 1,domx(2), domy(2), aux3
-    integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init, cpu_countrate,ic1, force_lingradT
+    integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init, cpu_countrate,ic1, force_lingradT, aux4, domx_aux(2)
     integer, target :: k
     real(dp), dimension(:,:), allocatable :: v, x, celula, rFUp !força n e n+1, velocidades e posições, [r*F, potencial e momento]
     real(dp), dimension(:), allocatable :: icell,jcell, nxv, nxv_send, nRfu, nRfu_send !dimensões das celulas e vetor de resultado pra imprimir
@@ -4142,6 +4156,8 @@ program main
     type(CFG_t) :: my_cfg
     character(5) :: particle
     character(4) :: wall
+    character(10) :: time 
+    character(8) :: date
     character(20) :: nome,arquivox, arquivov = '%' !nome de partícula deve ter até 20 caracteres
     type(string) :: part_nomes(10) ! vetor de strings
     character(1) :: optio
@@ -4260,6 +4276,11 @@ program main
     call CFG_get(my_cfg,"global%GField",GField)
     call CFG_get(my_cfg,"global%print_TC",print_TC)
 
+    if (termostato_vel_scaling .and. force_lingradT > 1) then
+        temp_Td(2) = t_fim
+        print*, "O termostato estará ligado ao longo de toda simulação"
+        print*, "temp_Td = ", temp_Td
+    end if
 
 
     if (NMPT < 1) NMPT = N
@@ -4527,10 +4548,14 @@ program main
 
     ! id = 0 é o processo raiz 
 
-    ! if (id == 0) then
+    if (id == 0) then
+        call date_and_time(date,time)
+        write(*,'("|| Program started at: ", a,":",a,":",a,2x, a,"/",a,"/",a, "||")') & 
+            time(1:2),time(3:4),time(5:8),date(5:6),date(7:8),date(1:4)
+
     !     print*,'Press Return to continue'
     !    read(*,*)
-    ! end if
+    end if
 
     if (id == 0) then 
         call system_clock(ic1,cpu_countrate)
@@ -4805,12 +4830,17 @@ program main
             end if 
 
             if (force_lingradT > 0 .and. i == interv_Td(j2)) then
-                do aux4 = 0, force_lingradT/subx
-
-                    domx_aux = [ (domx(2)-domx(1))/(force_lingradT/subx)
-
-                    call call corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,t)
-
+                ! ao longo da direção x, o gradiente terá [force_lingradT] passos. Teremos (domx(2)-domx(1))/(force_lingradT/subx) celulas por passo
+                do aux4 = 0, force_lingradT/subx-1
+                                      !cellx por processo  ! passos por processo 
+                    domx_aux = [1+aux4*(domx(2)-domx(1)) / (force_lingradT/subx), & 
+                                (aux4+1)*(domx(2)-domx(1)) / (force_lingradT/subx) ]  
+                    if (domx_aux(1) == 1) domx_aux(1) = 2
+                    if (domx_aux(2) == mesh(1)+2) domx_aux(2) = mesh(1) + 1 
+                    Td = sum(domx_aux)*(Td_hot-Td_cold)/(2*mesh(2))
+                    call corr_Kglobal(malha,domx_aux,domy,Td,propriedade, np,id,t,.true.)
+                    ! Este .true. indica para ignorar mpi, ou seja, pedaços de outras regiões 
+                    ! não serão levados em conta.
 
                     j2 = j2 + 1
                 end do
@@ -4980,7 +5010,9 @@ program main
         finish = real(ic1)/real(cpu_countrate,kind(0.d0))
         print '("Time = ",f10.3," seconds.")',(finish-start)
         open(unit=22,file='settings.txt',status="old", position="append", action="write")
-        write(22,*) "#:Execution time = ",(finish-start)," seconds."
+        write(22,'("#:Time to complete: ",f10.3," secounds.")') (finish-start)
+        write(22,'("#:Execution date: ",a,"/",a,"/",a,2x,a,":",a,":",a)') & 
+        date(5:6),date(7:8),date(1:4),time(1:2),time(3:4),time(5:8)
         close(22)
         call system('python csv2vtk_particles.py')
     end if

saida.f90

@@ -61,13 +61,13 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start)
             ! call cpu_time(time)    
                 call system_clock(ic1,cpu_countrate)
                 time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0))
-            if (step == 1) then 
+            if (step == 2) then 
                 call system_clock(ic1,cpu_countrate)
                 time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0))
                 etc = ((time - start)/real(step,kind(0.d0)) + (time-timep))*0.5*real(nimpre,kind(0.d0)) - (time - start)
                 dtimepp = (time-timep)
                 print '("Salvo arquivo ", A, "  t = ", f10.3, "  ETC: ", f10.3, "s" )',prop//extensao//'.'//trim(passo),t,etc                                    
-            else if (step > 1) then
+            else if (step > 2) then
                 call system_clock(ic1,cpu_countrate)
                 time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0))
                 etc = ((time - start)*6/real(step,kind(0.d0)) + ((time-timep)*2 + dtimepp*2))* & 

settings.ini

@@ -4,26 +4,26 @@
 # Argonio: epsilon/kb = 120K, sigma = 0.341nm 
 [global]
     nimpre  = 100 # quntidade de saidas 
-    N = 584 # número de partículas 
+    N = 5100  # número de partículas 
     Ntype = 2 # número de tipos de partícula
-    t_fim = 50  # tempo de execução 
+    t_fim = 10000  # tempo de execução 
     nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). 
-    dt = 0.0001 # passo de tempo
-    dimX = 400 # Dimensões da região de cálculo
-    dimY = 400 #
-    mesh = 80 80 # 300 60 # malha(elementos)
+    dt = 0.001 # passo de tempo
+    dimX = 20000 # Dimensões da região de cálculo
+    dimY = 5000 #
+    mesh = 200 40 # 300 60 # malha(elementos)
     rcut = 3 # *sigma
     wall = 'ppee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste
     termostato_nose_hoover = .false. #liga termostato Nose Hoover. Este funciona durante toda execução.
-    termostato_vel_scaling = .false.  # liga termostato por velocity scaling
+    termostato_vel_scaling = .true.  # liga termostato por velocity scaling
     Mc = 1
     Td = -1 #temperatura desejada global constante (-1 = não aplicar) kb = 1.38064852E-23, Td = (2/(3*kb))*E/epsilon(literatura)
-    temp_Td = 0 10 100 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] 
-    cold_cells = 1 10 1 40  # [cellx_ini, cellx_fin, celly_ini, celly_fin] banho frio
-    hot_cells = 190 200 1 40 # [cellx_ini, cellx_fin, celly_ini, celly_fin] banho quente
-    Td_hot = 1 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato.
-    Td_cold = 0.2 # v = sqrt(2*T) temperatura nas hot cells constante. Td_cold = -1 desliga termostato
-    force_lingradT = 10 # se > 0 irá dividir região no numero de subregiões entradas na direação X e será forçado um gradiente de temperatura linear. A região dos banhos será ignorada. A temperatura será controlado por um termostato velocity scaling
+    temp_Td = 0 35000 250 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] 
+    cold_cells = 1 10 1 40  # [cellx_ini, cellx_fin, celly_ini, celly_fin]
+    hot_cells = 190 200 1 40 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
+    Td_hot = 2 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato
+    Td_cold = 1 # v = sqrt(2*T) temperatura nas hot cells constante. Td_cold = -1 desliga termostato
+    force_lingradT = 20 # (NUMERO MULTIPLO DE 4 RECOMENDADO!) se > 0 irá dividir região no numero de subregiões entradas na direação X e será forçado um gradiente de temperatura linear. A região dos banhos será ignorada. A temperatura será controlado por um termostato velocity scaling
     vd = 0 0 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') 
     NMPT = 1000 # Número máximo de partículas que poderá mudar de process. Se não souber estimar, usar -1 ou o valor de N.
     GField = 0 0 # Campo de força gravitacional que afeta todas as partículas
@@ -33,29 +33,29 @@
 # se v_file é dada, então ela será usada para dar a velocidade às partículas. Deve ser arquivo .csv. Deixar "%%" se não usar. 
 # se ismolecule = .true. então a partícula é levada em conta no termostato 
 [par_0]
-    x = 'poucas1.csv' 
+    x = 'moleculas_sc.csv' 
     v = 0 0 # velocidade global
     v_file = 'v_file_0.csv' # velocidade de cada partícula
     m = 1
     nome = 'g1'
     sigma = 1 # Lennard Jones
     epsilon = 1 # Lennard Jones
-    quantidade = 484
+    quantidade = 5000
     x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay
     rs = 0 # raio sólido. Posição de partículas na superfície de um sólido de raio rs
     fric_term = 0 # fricção artifical não funciona quando o termostato nose hoover está ligado
     ismolecule = .true.			
 [par_1]
-    x = 'poucas2.csv' 
+    x = 'particulas_sc.csv' 
     v = 0 0 # velocidade global
     v_file = '%%v_file_0.csv' # velocidade de cada partícula
-    m = 50
+    m = 100
     nome = 'g2'
-    sigma = 1 # Lennard Jones
-    epsilon = 1 # Lennard Jones
+    sigma = 0.94 # Lennard Jones
+    epsilon = 0.428 # Lennard Jones
     quantidade = 100
-    x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay
-    rs = 1 # raio sólido. Posição de partículas na superfície de um sólido de raio rs
+    x_lockdelay = 1000 # só vai poder mudar de posição a partir de t = x_lockdelay
+    rs = 10 # raio sólido. Posição de partículas na superfície de um sólido de raio rs
     fric_term = 0 # fricção artifical	
     ismolecule = .false.						
 #[par_2]