ignorar particulas em compK implementado

csv2vtk_particles.py

@@ -57,7 +57,10 @@
 N = int(config['global']['N'].split()[0])
 nimpre =  int(config['global']['nimpre'].split()[0])
 ntype = int(config['global']['Ntype'].split()[0])
-print_TC = int(config['global']['print_TC'].split()[0])
+if config['global']['print_TC'].split()[0] == '.false.':
+    print_TC = False
+else:
+    print_TC = True
 quant = []
 rs = [] # raio sólido
 sigma = []
@@ -107,7 +110,7 @@
 
 print('Converting...')
 if pbar:
-    bar = progressbar.ProgressBar(max_value=nimpre-1)
+    bar = progressbar.ProgressBar(max_value=nimpre)
 for fnum in range(0,nimpre+1):
     with open('temp/position.csv.'+str(fnum),encoding='utf-8') as file_locus:
         csv_lector = csv.reader(file_locus,delimiter = ',')
@@ -156,7 +159,7 @@
 
 shutil.rmtree('temp')
 
-if print_TC == 1:
+if print_TC:
     a = os.listdir('temp2')
     zip_rFup = ZipFile(via+'/'+folder+'/rFuP.zip','w')
     for f in a:

data.f90

@@ -13,7 +13,6 @@ module data
         real(dp),dimension(2) :: F !força nela
         integer               :: grupo !grupo que a partícula pertence
         integer               :: n !numero da partícula
-        ! integer, dimension(2) :: mic !minimal image convenson (ajuda a contar quantas vezes a partcula atravessou a borda numa direção)
         logical               :: flag ! bandeira auxiliar
     end type data_t
 

debug_lennard.sh

@@ -0,0 +1,2 @@
+#!/bin/bash
+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

lennard.f90

@@ -491,7 +491,7 @@ function comp_Kglobal(malha,domx,domy,propriedade,np,id,t) result(K)
                 do while (associated(node))
                     ptr = transfer(list_get(node), ptr)
                     !calcula a energia cinética atual
-                    if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then
+                    if (propriedade(ptr%p%grupo)%x_lockdelay <= t  .and. propriedade(ptr%p%grupo)%ismolecule) then
                         K = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + K
                         nump = nump +1
                     end if
@@ -544,7 +544,7 @@ subroutine corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,tt)
                 node => list_next(malha(i,j)%list)
                 do while (associated(node))
                     ptr = transfer(list_get(node), ptr)
-                    if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then
+                    if (propriedade(ptr%p%grupo)%x_lockdelay <= t  .and. propriedade(ptr%p%grupo)%ismolecule) then
                         !calcula a energia cinética atual
                         ptr%p%v = beta*ptr%p%v ! aqui o T é o beta 
                     end if
@@ -579,7 +579,7 @@ function comp_K(malha,domx,domy,s_cells,propriedade,np,id,t) result(K)
                     do while (associated(node))
                         ptr = transfer(list_get(node), ptr)
                         !calcula a energia cinética atual
-                        if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then
+                        if (propriedade(ptr%p%grupo)%x_lockdelay <= t .and. propriedade(ptr%p%grupo)%ismolecule) then
                             K = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + K
                             nump = nump +1
                         end if
@@ -636,7 +636,7 @@ subroutine corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedad
                     node => list_next(malha(i,j)%list)
                     do while (associated(node))
                         ptr = transfer(list_get(node), ptr)
-                        if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then
+                        if (propriedade(ptr%p%grupo)%x_lockdelay <= t  .and. propriedade(ptr%p%grupo)%ismolecule) then
                             !calcula a energia cinética atual
                             ptr%p%v = betac*ptr%p%v ! aqui o T é o beta 
                         end if
@@ -2546,7 +2546,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, print_TC
+    integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init
     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
@@ -2564,7 +2564,7 @@ program main
     type(string) :: part_nomes(10) ! vetor de strings
     character(1) :: optio
     type(data_ptr) :: ptr !integer,pointer :: ptr,ptrn
-    logical :: laux, termostato
+    logical :: laux, termostato, print_TC
     character(len=32) :: arg
     character :: arg1(1)
     type(lstr) :: LT
@@ -2637,7 +2637,7 @@ program main
         "Max Number of particles that will change process per iteraction")
     call CFG_add(my_cfg,"global%GField",(/0.0_dp, 0.0_dp/), &
         "Uniform Gravitational Field")
-    call CFG_add(my_cfg,"global%print_TC",1, &
+    call CFG_add(my_cfg,"global%print_TC",.false., &
         "Print transport coefficient data")
 
     call CFG_read_file(my_cfg, "settings.ini")
@@ -2711,13 +2711,14 @@ program main
         call CFG_add(my_cfg, particle//"%rs",1.1_dp, &
             "solid radius "//particle) 
         call CFG_add(my_cfg,particle//"%fric_term",1.0_dp, &
-            "Friction term")    
+            "Friction term")  
+        call CFG_add(my_cfg,particle//"%ismolecule",.true., &
+            "Is Molecule")      
     end do
 
     dx_max = 10*dimx !pra definir critério de estabilidade no uso de malha
     do i = 0,(Ntype-1)
         write(particle,'(a,i0)') 'par_',i
-
         call CFG_get(my_cfg, particle//"%quantidade", quant)
         call CFG_get(my_cfg, particle//"%nome", nome)
         part_nomes(i+1)%str = nome
@@ -2729,7 +2730,8 @@ program main
         call CFG_get(my_cfg, particle//"%x_lockdelay", propriedade(i+1)%x_lockdelay)
         call CFG_get(my_cfg, particle//"%rs", propriedade(i+1)%rs)
         call CFG_get(my_cfg, particle//"%fric_term",propriedade(i+1)%fric_term)
-
+        call CFG_get(my_cfg, particle//"%ismolecule",propriedade(i+1)%ismolecule)
+
         ! le o arquivo com posições
         if (dx_max < propriedade(i+1)%sigma*rcut/2) then
             dx_max = propriedade(i+1)%sigma*rcut/2
@@ -2965,7 +2967,7 @@ program main
         i = interv(j)
     end if
 
-    if (print_TC == 1) then
+    if (print_TC) then
         if (wall(1:2) == 'pp' .or. wall(3:4) == 'pp') then
             allocate(nRfu(N*5),nRfu_send(N*6),rFUp(N,7))
         else 
@@ -3060,7 +3062,7 @@ program main
 
             deallocate(nxv,nxv_send)
 
-            if (print_TC == 1) then 
+            if (print_TC) then 
                 if ((wall(1:2) == 'pp' .or. wall(3:4) == 'pp') .and. id /= 0) then
                     mic_trf(1:N) = mic(:,1)
                     mic_trf(N+1:2*N) = mic(:,2)

mod0.f90

@@ -14,6 +14,7 @@ module mod0
             real(dp) :: rs !raio sólido 
             real(dp)  :: x_lockdelay    ! só vai poder mudar de posição a partir de t = x_lockdelay
             real(dp) :: fric_term
+            logical :: ismolecule
         end type prop_grupo
 
             ! LSTR  lista de transferência pro MPI

pack.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+read -e -p "Enter particle generator file: " VAR1
+
+if [ "$VAR1" = "" ]; then
+    echo "No particle generator given."
+    tar -cvf lennard_files.tar csv2vtk_particles.py lennard settings.ini verify_settings.py visualizar.py
+else
+    tar -cvf lennard_files.tar csv2vtk_particles.py lennard $VAR1 settings.ini verify_settings.py visualizar.py
+fi
+

settings.ini

@@ -23,22 +23,24 @@
     vd = 2 2 # 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
-    print_TC = 1 # imprimir dados para calcular coeficiente de transporte
+    print_TC = .false. # imprimir dados para calcular coeficiente de transporte
 # PARTICLE PROPERTIES
 # v é velocidade global, entrar dois número correspondentes ao vetor vx vy que serão aplicadas em todas as partículas
-# 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 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 = 'gas.csv' 
     v = 0 0 # velocidade global
-    v_file = '%v_file_0.csv' # velocidade de cada partícula
+    v_file = '%%v_file_0.csv' # velocidade de cada partícula
     m = 1
     nome = 'g1'
     sigma = 1 # Lennard Jones
     epsilon = 1 # Lennard Jones
     quantidade = 2000
     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				
+    fric_term = 0 # fricção artifical	
+    ismolecule = .true.			
 #[par_1]
 #    x = 'bin_peq2_5x1.csv' 
 #    v = 0 0 # velocidade global
@@ -62,4 +64,4 @@
 #    quantidade = 20
 #    x_lockdelay = 300 # só vai poder mudar de posição a partir de t = x_lockdelay
 #    rs = 5
-#    fric_term = 0 # fricção artifical		
+#    fric_term = 0 # fricção artifical		

verify_settings.py

@@ -1,5 +1,6 @@
 
 import configparser
+import os
 
 change = False
 aux = True
@@ -26,9 +27,12 @@
     try:
         quant.append(int(config['par_'+str(i)]['quantidade'].split()[0])) 
         print("Number of particle types greater then expected (ntype = {} but there is par_0 to par_{}).\n".format(ntype,i))
-        nty = input("Enter the right number of particle types.\n")
-        i = int(nty) 
-        config.set('global','Ntype',nty)
+        os.system('xdg-open settings.ini')
+        input("Fix the config file, then enter to continue...")
+        config.read('settings.ini')
+        # nty = input("Enter the right number of particle types.\n")
+        # i = int(nty) 
+        # config.set('global','Ntype',nty)
     except KeyError:
         aux = False 
         pass 
@@ -41,44 +45,52 @@
     try:
         for i in range(ntype):
             quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))
-            j+= 1
         aux = False 
     except:
-        print("Number of particle types (ntype) smaller then expected. It will correct to proceed with ntype = {}".format(j))
-        config.set('global','Ntype',str(j))
-        ntype = j 
-        j = 0
-        change = True
+        print("Number of particle types (ntype) smaller then expected. It maybe correct to proceed with ntype = {}".format(j))
+        os.system('xdg-open settings.ini')
+        input("Fix the config file, then enter to continue...")
+        config.read('settings.ini')
+
+aux = True
+while aux:
+    for i in range(ntype):
+        x_file = config['par_'+str(i)]['x'].split()[0]
+        x_file = x_file.replace("'","")
+        x_file = x_file.replace('"','')
+        if os.path.isfile(x_file):
+            aux = False
+        else:
+            print("Position file {} not found.".format(x_file))
+            os.system('xdg-open settings.ini')
+            input("Fix the config file, then enter to continue...")
+            config.read('settings.ini')
+        v_file = config['par_'+str(i)]['v_file'].split()[0]
+        v_file = v_file.replace("'","")
+        v_file = v_file.replace('"','')
+        if v_file[0] == '%' or os.path.isfile(x_file):
+            aux = False
+        else:
+            print("Velocity file {} not found.".format(v_file))
+            os.system('xdg-open settings.ini')
+            input("Fix the config file, then enter to continue...")
+            config.read('settings.ini')
 
 
 if sum(quant) != N: 
-    print("The particle groups have a different number of particles then specified {}. \nIt will be corrected to {}".format(N,sum(quant)))
-    config.set('global','N',str(sum(quant)))
-    change = True
-
-if change:
-    print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n")
-    print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n")
-    a = input("Changes will be made to the settings file. Enter to proceed or Ctrl+C to abort.")
-    a = input()
-    with open('settings.ini', 'w') as configfile:
-        config.write(configfile)
+    print("The particle groups have a different number of particles then specified {}. \nIt should maybe be corrected to {}".format(N,sum(quant)))
+    os.system('xdg-open settings.ini')
+    input("Fix the config file, then enter to continue...")
+    config.read('settings.ini')
+
 
 try:
     from evtk.hl import pointsToVTK 
 except ModuleNotFoundError:
     print("Evtk module not found! You will not be able to convert the .CSV files to .VTU. Try:\n")
     print("conda install -c e3sm evtk     OR\npip install pyevtk")
 
-with open('settings.ini','r') as file:
-    data = file.readlines()
-
-for line in range(len(data)):
-    if data[line][0] != "[" and data[line][0] != "#" and data[line][0] != " ":
-        data[line] = '\t' + data[line]
 
-with open('settings.ini','w') as file:
-    file.writelines(data)