diff --git a/csv2vtk_particles.py b/csv2vtk_particles.py new file mode 100644 index 0000000..a29a57e --- /dev/null +++ b/csv2vtk_particles.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jan 9 13:03:06 2019 + +Este programa lê arquivos CSV onde cada arquivo é um passo de tempo e cada +linha é uma partícula e exporta em vtk o reultado. + +@author: gabriel +""" + +from evtk.hl import pointsToVTK +import numpy as np +import csv, os, shutil +import configparser + +via = os.getcwd() + +aux = True +folder = 'result' +i = 1 + +while aux: + try: + os.mkdir(folder) + aux = False + except FileExistsError: + print('Folder {} exists'.format(folder)) + opt = input('Overwrite? [y/n]') + if opt == 'y' or 'opt' == 'Y': + aux =False + else: + folder = input('Enter new folder name:\n') + if len(folder) == 0: + folder = 'result' + '_' + str(i) + aux = True + i = i+1 + +print('The results will be saved at {}'.format(folder)) + +#Ler a quantidade de arquivos de settings. + +config = configparser.ConfigParser() +config.read('settings.ini') + +N = int(config['global']['N'].split()[0]) +nimpre = int(config['global']['nimpre'].split()[0]) +ntype = int(config['global']['Ntype'].split()[0]) +quant = [] + +for i in range(ntype): + quant.append(int(config['par_'+str(i)]['quantidade'].split()[0])) + +a = os.listdir('temp') +if len(a)/2 < nimpre: + print("Propable incomplete execution. Processing {} files\n.".format(len(a)/2)) + nimpre = int(len(a)/2)-1 + +tipo = np.zeros(N) +j,k = 0,0 +for i in range(len(quant)): + for j in range(quant[i]): + tipo[j+k] = i + k = quant[i] + +# Ler os arquivos e colocá-los em vetores + +x = np.zeros(N) +y = np.zeros(N) +z = np.zeros(N) + +vx = np.zeros(N) +vy = np.zeros(N) +vz = np.zeros(N) + +cx = np.zeros(N) +cy = np.zeros(N) +cz = np.zeros(N) + +try: + shutil.move(via+'/settings.txt',via+'/'+folder+'/settings.txt') +except: + print('No settings file found!\n') + +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 = ',') + i = 0 + for linea in csv_lector: + x[i] = linea[0] + y[i] = linea[1] + i = i+1 + shutil.move(via+'/temp/position.csv.'+str(fnum),via+'/'+folder+'/position.csv.'+str(fnum)) + + with open('temp/velocity.csv.'+str(fnum),encoding='utf-8') as file_velocitas: + csv_lector = csv.reader(file_velocitas,delimiter = ',') + i = 0 + for linea in csv_lector: + vx[i] = linea[0] + vy[i] = linea[1] + + i = i+1 + + shutil.move(via+'/temp/velocity.csv.'+str(fnum),via+'/'+folder+'/velocity.csv.'+str(fnum)) + + #with open('cell.csv.'+str(fnum),encoding='utf-8') as file_velocitas: + #csv_lector = csv.reader(file_velocitas,delimiter = ',') + #i = 0 + #for linea in csv_lector: + #cx[i] = int(float(linea[0])) + #cy[i] = int(float(linea[1])) + + #i = i+1 + + #shutil.move(via+'/cell.csv.'+str(fnum),via+'/'+folder+'/cell.csv.'+str(fnum)) + + pointsToVTK(via+'/'+folder +'/result_'+str(fnum), x, y, z, data = {"Vx" : vx, "Vy" : vy, "Tipo" : tipo}) #, "Cellx" : cx, "Celly" : cy + +shutil.rmtree('temp') + + diff --git a/data.f90 b/data.f90 new file mode 100644 index 0000000..7224afa --- /dev/null +++ b/data.f90 @@ -0,0 +1,31 @@ +module data + use mod1 + use linkedlist + + private + public :: data_t + public :: data_ptr + + ! Data is stored in data_t + type :: data_t + real(dp),dimension(2) :: x !posição + real(dp),dimension(2) :: v !velocidade + real(dp),dimension(2) :: F !força nela + integer :: grupo !grupo que a partícula pertence + integer :: n !numero da partícula + logical :: flag ! bandeira auxiliar + end type data_t + + ! A trick to allow us to store pointers in the list + type :: data_ptr + type(data_t), pointer :: p + end type data_ptr + + + + ! pra fazer vetor de ponteiro +! type :: container +! type(list_t), pointer :: list => null() +! end type container + +end module data diff --git a/lennard.f90 b/lennard.f90 new file mode 100644 index 0000000..aa9a758 --- /dev/null +++ b/lennard.f90 @@ -0,0 +1,2360 @@ + +module mod0 + use mod1 + use linkedlist + type :: container + type(list_t), pointer :: list => null() + end type container + + ! Propriedades de cada grupo de partícula (indicado pelo índice do vetor com este tipo) + ! nos permite mais fases e economiza espaço + type :: prop_grupo + real(dp) :: m ! mass + real(dp) :: epsilon + real(dp) :: sigma + real(dp) :: x_lockdelay ! só vai poder mudar de posição a partir de t = x_lockdelay + end type prop_grupo + + ! LSTR lista de transferência pro MPI + type :: lstr + real(dp),allocatable :: lstrdb_N(:) + real(dp),allocatable :: lstrdb_S(:) + real(dp),allocatable :: lstrdb_E(:) + real(dp),allocatable :: lstrdb_W(:) + real(dp),allocatable :: lstrdb_D(:) ! Diagonal. Para transferir nas direções NE, NS, SE, SW no encontro de 4 processos + + integer,allocatable :: lstrint_N(:) + integer,allocatable :: lstrint_S(:) + integer,allocatable :: lstrint_E(:) + integer,allocatable :: lstrint_W(:) + integer,allocatable :: lstrint_D(:) + end type lstr + + +end module mod0 + +module mod2 + !use mod0 + !use data + contains + + !Calcula velocidades distribuidas de acordo com MaxwellBolzmann + ! factor = sqrt() = sqrt(kb*T/m) por componente + subroutine MaxwellBoltzmann(malha,mesh,factor) + use linkedlist + use mod1 + use randnormal + use mod0 + use data + + real(dp),intent(in) :: factor(2) + integer :: i,j + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer, intent(in) :: mesh(2) + type(list_t), pointer :: node + type(data_ptr) :: ptr + + do i = 2,mesh(2)+1 + do j = 2,mesh(1)+1 + ! print *, 'posição', i, ',', j + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + ptr%p%v(1) = ptr%p%v(1) + factor(1)*GaussDeviate() + ptr%p%v(2) = ptr%p%v(2) + factor(2)*GaussDeviate() + ptr%p%F = [0,0] + node => list_next(node) + end do + end do + end do + + end subroutine MaxwellBoltzmann + + !computa o potencial total + function comp_pot(mesh,malha,propriedade,r_cut) result(V) + use linkedlist + use mod1 + use data + use mod0 + type(data_ptr) :: ptr, ptrn + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + real(dp) :: rcut,r_cut + integer :: i,j + integer, intent(in) :: mesh(:) + real(dp) :: r, aux2(2),V, x1(2),x2(2),sigma, epsil + type(list_t), pointer :: node, next_node + V = 0 + do i = 1,mesh(1)+2 + do j = 1,mesh(2)+2 +! ! !print*, 'Linha 67',i,j + ! dentro da malha(i,j) + node => list_next(malha(i,j)%list) +! print*,'B',i,j +! print*,associated(node),i,j + do while (associated(node)) + ptr = transfer(list_get(node), ptr) !particula selecionada + x1 = ptr%p%x +! ! !print*, 'L73' +! !read(*,*) + !calcular a força desta com todas as outras partículas + next_node => list_next(node) ! próxima partícula da célula +! ! !print*, 'L77',ptr + node => list_next(node) ! a ser computado com a particula selecionada + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + x2 = ptrn%p%x + sigma = propriedade(ptr%p%grupo)%sigma + epsil = propriedade(ptr%p%grupo)%epsilon + rcut = sigma*r_cut + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + V = abs((sigma/r)**6*((sigma/r)**6-1)) + V + end if + node => list_next(node) ! próxima partícula da célula + end do + + !Células ao redor + if (i /= mesh(1)+2) then + if (j == mesh(2)+2) then + node => list_next(malha(i+1,1)%list) + else + node => list_next(malha(i,j+1)%list) + end if + + do while (associated(node)) + + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + x2 = ptrn%p%x + ! ! !print*, 'L120' + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + ! sigma = propriedade(ptr%p%grupo)%sigma + ! epsilon = propriedade(ptr%p%grupo)%epsilon + ! rcut = sigma*rcut + if (r <= rcut) then + V = abs((sigma/r)**6*((sigma/r)**6-1)) + V + end if + node => list_next(node) ! próxima partícula da célula + + end do + + end if + + end do + ! print*,'linha137' + end do +! print*,'L143 LJ' + end do + + V = V*4*epsil + + end function comp_pot + + function comp_K(malha,domx,domy,propriedade,np,id,t) result(K) + use linkedlist + use mod1 + use data + use mod0 + use mpi + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + integer, intent(in) :: domx(2),domy(2) + type(container), allocatable,dimension(:,:),intent(in) :: malha + type(data_ptr) :: ptr + real(dp) :: kb = 1.38064852E-23,K, auxres(8), nump = 0, aux(2) + type(list_t), pointer :: node + integer ( kind = 4 ) :: np, ierr, id + real(dp), intent(in) :: t + + do i = domy(1),domy(2) + do j = domx(1),domx(2) + ! print *, 'posição', i, ',', j + node => list_next(malha(i,j)%list) + 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 + K = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + K + nump = nump +1 + end if + node => list_next(node) + end do + end do + end do + + ! então o processo é paralelo, vamos juntar tudo + aux = [K, nump] + if (np > 1) then + call MPI_GATHER(aux2, 2, MPI_DOUBLE_PRECISION, Kres, 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (id == 0) then + K = auxres(1)+auxres(3)+auxres(5)+auxres(7) + nump = auxres(2)+auxres(4)+auxres(6)+auxres(8) + auxres = [K,nump,K,nump,K,nump,K,nump] + end if + call MPI_SCATTER(auxres, 2, MPI_DOUBLE_PRECISION, aux, 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + end if + K = aux(1)/aux(2) + + end function comp_K + + !Corrige energia cinética/configura temepratura + subroutine corr_K(malha,domx,domy,N,Td,propriedade, np,id,tt) + use linkedlist + use mod1 + use data + use mod0 + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + real(dp),intent(in) :: Td!energia cinética + real(dp) :: T,K,kb = 1.38064852E-23 + integer, intent(in) :: N,domx(2),domy(2) + type(container), allocatable,dimension(:,:),intent(in) :: malha + type(data_ptr) :: ptr + type(list_t), pointer :: node + integer ( kind = 4 ) :: np, id + real(dp), intent(in) :: tt + !calcula temperatura atual + T = (2/(3*kb))*comp_K(malha,domx,domy,propriedade,np,id,tt) + T = sqrt(Td/T) ! aqui o T é o Beta + do i = domy(1),domy(2) + do j = domx(1),domx(2) + ! print *, 'posição', i, ',', j + node => list_next(malha(i,j)%list) + ! print*, i,j + ! print*, 'is associated?', associated(node) + ! print*, associated(node) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + !calcula a energia cinética atual + ptr%p%v = T*ptr%p%v ! aqui o T é o beta + end if + node => list_next(node) + end do + end do + end do + end subroutine corr_K + + !força de ficção + function comp_fric(r,v,fric_term) result(f) + use mod1 + real(dp), intent(in) :: r(2),v(2),fric_term + real(dp) :: f(2) + + f = (fric_term/(r(1)**2 + r(2)**2))*(v(1)*r(1)+v(2)*r(2))*[r(1),r(2)] + end function comp_fric + + + ! atualiza forças + subroutine comp_F(GField, mesh,malha,propriedade,r_cut,fric_term,domx,domy, ids, id, t) + use linkedlist + use mod1 + use data + use mod0 + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + real(dp), intent(in) :: GField(2), t + real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção + real(dp) :: x1(2),v1(2),x2(2),v2(2) + integer :: i,j,ct = 0 !,ptr, ptrn + integer, intent(in) :: mesh(:),domx(2),domy(2) + real(dp) :: Fi(2)=0,r, aux2(2),fR(2) + type(list_t), pointer :: node, next_node + type(data_ptr) :: ptr,ptrn + integer ( kind = 4 ), intent(in) :: ids(8) + + + !Lennard Jones + fR = [0,0] + do i = domy(1),domy(2) ! i é linha + do j = domx(1),domx(2) + node => list_next(malha(i,j)%list) + if (associated(node)) then + ptr = transfer(list_get(node), ptr) + ! print*, "L 253", ptr%p%n + end if + do while (associated(node)) + ! print*, "Encontrada", ct, "celulas, id=", id + ptr = transfer(list_get(node), ptr) !particula selecionada + ptr%p%flag = .true. ! indica ao comp_x que a partícula precisa ser calculada + x1 = ptr%p%x + v1 = ptr%p%v + m1 = propriedade(ptr%p%grupo)%m + sigma_a = propriedade(ptr%p%grupo)%sigma + ! rcut = r_cut*sigma + epsil_a = propriedade(ptr%p%grupo)%epsilon + if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + ptr%p%F = ptr%p%F + GField*m1 + end if + ! print '("x1 =", f10.6, " ", f10.6, " n ", i2, " cell ",i2," ",i2)',x1(1),x1(2), ptr%p%n,i,j + ! if (id == 0) read(*,*) + !calcular a força desta com todas as outras partículas + next_node => list_next(node) ! próxima partícula da célula + node => list_next(node) ! a ser computado com a particula selecionada + + ! NA PRÓPRIA CELULA + do while (associated(node)) + ! print*, "x2", x2 + ! print*, "ID", id + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a *epsil_b ) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a *epsil_b ) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + ! print '("r ", f10.5, " rcut", f10.5)',r,rcut + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + !! PARA CASO PARALELO, INTERAGIR COM AS CELULAS AO NORTE E A OESTE + ! + ! SUL + if (ids(2) > -1 .and. i == domy(1)) then + node => list_next(malha(i-1,j)%list) !interagirá com a anterior linha apenas + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ! ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i-1,j-1)%list) !interagá com a linha e coluna anterior + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ! ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + + if (j /= mesh(1)+2) then !se não for a última coluna + node => list_next(malha(i-1,j+1)%list) !interagá com a linha e coluna anterior + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ! ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + end if + + ! OESTE + if (ids(4) > -1 .and. j == domx(1)) then + node => list_next(malha(i,j-1)%list) !interagirá com a anterior linha apenas + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ! ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (i /= mesh(2)+2) then !se não for a primeira linha + node => list_next(malha(i+1,j-1)%list) !interagá com a linha e coluna anterior + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ! ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + end if + !! fim do caso para paralelismo + + !Células ao redor !i é linha, j é coluna + if (i /= mesh(2)+2) then ! se não for a última linha + if (j == mesh(1)+2) then ! se for a última coluna + node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha apenas + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + else + !interagirá com a próxima linha e coluna, e na diagonal + node => list_next(malha(i,j+1)%list) + + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + node => list_next(malha(i+1,j+1)%list) !interagirá com a próxima linha e coluna + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil_a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + end if + + else ! se for a última lina, só interage com a celua ao lado + node => list_next(malha(i,j+1)%list) !interagirá com a próxima linha e coluna + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + epsil_b = propriedade(ptrn%p%grupo)%epsilon + ! Lorenz-Betherlot rule for mixing epsilon sigma + if (sigma_a > sigma_b) then + rcut = r_cut*sigma_a + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + sigma = sigma_a + epsil = epsil__a + end if + + x2 = ptrn%p%x + v2 = ptrn%p%v + r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) !raio + if (r <= rcut) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*[(x1(1)-x2(1)), (x1(2)-x2(2))] + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1)),-(x1(2)-x2(2))],[-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + ! ! ! print*, "L F431", ptr%p%F, id,i,j + node => next_node + end do + ! print*,'linha137' + end do + ! print*,'L143 LJ' + end do + ! print*, "x_1_LJ", x1 + !print*, 'Linha 143' + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end subroutine comp_F + + ! atualiza posições + subroutine comp_x(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,domy,id, np) + use linkedlist + use mod1 + use data + use mod0 + use mpi + use matprint + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer :: i,j,k, cell(2), status(MPI_STATUS_SIZE), count, destD + integer, intent(in) :: N,mesh(2),domx(2),domy(2) + type(list_t), pointer :: node, previous_node + real(dp), intent(in) :: dt, t + type(data_ptr) :: ptr + real(dp) :: x(2),m + real(dp), intent(in) :: icell(:), jcell(:) + integer ( kind = 4 ), intent(in) :: ids(8), id, np + integer ( kind = 4 ) :: tag,cont_db(8),cont_int(8) + type(lstr) :: LT + ! real(dp),intent(inout) :: celula(:,:) + ! IDS correspondem às celulas [N,S,E,W] + ! !print*, 'L SOMA', (ids(1) + ids(2) + ids(3) +ids(4)), 'ID', id + !!! ! print*, "L 464", id + cont_db = 0 + cont_int = 0 + ! if (id == 0) read(*,*) + + ! Esvazia as celulas emprestadas + if (domy(1) > 1) then + i = domy(1)-1 + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "A", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + ! node => list_next(malha(i,j)%list) + ! if (associated(node)) then + ! ! print*, "LIMPANDO A", i,j + ! ! ! print*, "L AA", ptr%p%n + ! ! print*, "AA", id + ! call list_free(malha(i,j)%list) + ! ! print*, "AAA", associated(node) + ! call list_init(malha(i,j)%list) + ! end if + end do + end if + if (domy(2) < mesh(2)+2) then + i = domy(2)+1 + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "B", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + ! node => list_next(malha(i,j)%list) + ! if (associated(node)) then + ! ! print*, "LIMPANDO B", i,j + ! ! print*, "BB", id + ! call list_free(malha(i,j)%list) + ! ! print*, "BBB", associated(node) + ! call list_init(malha(i,j)%list) + ! end if + end do + end if + if (domx(2) < mesh(1)+2) then + j = domx(2)+1 + do i = domy(1),domy(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "D", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + ! node => list_next(malha(i,j)%list) + ! if (associated(node)) then + ! ! print*, "LIMPANDO C", i,j + ! ! print*, "DD", id + ! call list_free(malha(i,j)%list) + ! ! print*, "DDD", associated(node) + ! call list_init(malha(i,j)%list) + ! end if + end do + end if + if (domx(1) > 1) then + j = domx(1)-1 + do i = domy(1),domy(2) + node => list_next(malha(i,j)%list) + previous_node => malha(i,j)%list + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + ! node => list_next(malha(i,j)%list) + ! if (associated(node)) then + ! ! print*, "LIMPANDO D", i,j + ! ! print*, "CC", id + ! ! ptr = transfer(list_get(node), ptr) + ! ! print*, ptr%p%n + ! call list_free(malha(i,j)%list) + ! ! print*, "CCC", associated(node) + ! call list_init(malha(i,j)%list) + ! end if + end do + end if + ! PARA 4 processos. Limpa o que está na diagonal adjacente + if (np > 1 .and. sum(ids(5:8)) > -4) then + if (domy(1) /= 1) i = domy(1) -1 + if (domy(2) /= mesh(2)+2) i = domy(2) +1 + if (domx(1) /= 1) j = domx(1) -1 + if (domx(2) /= mesh(1) + 2) j = domx(2) +1 + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + ! node => list_next(malha(i,j)%list) + ! print*, "Limpando malha", i,j, "ID", id + ! if (associated(node)) then + ! ! print*, "LIMPANDO E", i,j + ! ! print*, "EE", id + ! ! ptr = transfer(list_get(node), ptr) + ! ! print*, ptr%p%n + ! call list_free(malha(i,j)%list) + ! ! print*, "EEE", associated(node) + ! call list_init(malha(i,j)%list) + ! end if + end if + ! print*, "cont_int", cont_int + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + + do i = domy(1),domy(2) + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, i,j, "id", id + ! ptr%p%flag = true significa que ainda não foi computado para evitar loop + + if (associated(node)) then + ptr = transfer(list_get(node), ptr) + ! aqui vemos se já passou o tempo que as posições ficaram travadas. + if (propriedade(ptr%p%grupo)%x_lockdelay > t) then + ptr%p%flag = .false. + end if + end if + + do while (associated(node) .and. ptr%p%flag) + ! print*, "L PART", i,j, "id", id + m = propriedade(ptr%p%grupo)%m + ! print*, 'grupo =', ptr%p%grupo, 'massa=', m + ptr%p%flag = .false. + ! ptr = transfer(list_get(node), ptr) + ptr%p%x(1) = ptr%p%x(1) + dt*ptr%p%v(1) + ptr%p%F(1)*dt**2/(2*m) + ptr%p%x(2) = ptr%p%x(2) + dt*ptr%p%v(2) + ptr%p%F(2)*dt**2/(2*m) + ! print*,'F_0 = ',ptr%p%F(1),ptr%p%F(2), "id", id + ! print '("x_0 = [",f10.3,", ",f10.3, "] ", "i, j =", i2, " ", i2, " n ",i2 )',ptr%p%x(1),ptr%p%x(2),i,j, ptr%p%n!"id", id, "n", ptr%p%n + ! Ordena as celulas + x = ptr%p%x + + ! Teste se escapou para a esquerda + if (x(1) <= jcell(j-1)) then + cell = [i,j-1] + !para esquerda para cima + if (x(2) <= icell(i-1)) then + cell = [i-1,j-1] + !para esquerda e para baixo + else if (x(2) > icell(i)) then + cell = [i+1,j-1] + end if + + ! Teste se escapou para a direita + else if (x(1) > jcell(j)) then + cell = [i,j+1] + !para direita para cima + if (x(2) <= icell(i-1)) then + cell = [i-1,j+1] + !para direita e para baixo + else if (x(2) > icell(i)) then + cell = [i+1,j+1] + ! print*, "icell(i)", icell(i), "jcell(j)", jcell(j), "x=", x + ! ! print*, "L PUTAMERDA" + end if + + else if (x(2) <= icell(i-1)) then + ! print*,'!Teste se escapou para cima' + ! print*,x(2) ,'<', icell(i-1) + cell = [i-1, j] + !Teste se escapou para baixo + else if (x(2) > icell(i)) then + ! print*,'!Teste se escapou para baixo' + ! print*,x(2), '>', icell(i) + cell = [i+1,j] + else + cell = [i,j] + end if + + ! print*, "Contint(1)", cont_int(1) + + ! VERIFICAÇÃO SE MUDOUD DE DOMÍNIO + if (cell(1) /= i .or. cell(2) /= j) then !Mudou de celula + ! print*, "MUDOU" + ! if (id == 0) read(*,*) + ! Se a partícula chegar na fronteira do domínio que um processador + ! cuida, então esta partícula é colocada na lista linkada referente + ! a este célula com list_change. Além disso ela precisa ser + ! transferida para o outro processo vizinho para que ele possa + ! calcular sua influência no domínio. + ! Se a partícula ultrapassar o domínio, então ela é transferida + ! para o próximo processo e dealocada da lista com list_remove + ! ! ! print*, "L FORÇA", ptr%p%F, id + if (cell(2) >= domx(2) .and. cell(2) < mesh(1)+2 ) then + ! print*, "L 538", cell(1), cell(2),"part", ptr%p%n,"i,j", i, j + ! 6 elementos + LT%lstrdb_E(cont_db(3)+1:cont_db(3)+6) = [x(1),x(2), & + ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + ! 4 elementos + ! print*, "id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "id",id,"transferindo para o leste", LT%lstrint_E, "cont", cont_int(3) + ! if (cell(2) > domx(2)) then + ! ! print*, 'affs' + ! ! read(*,*) + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_E + + ! ! call list_remove(previous_node) + ! ! print*, 'affs' + ! ! node => list_next(previous_node) + ! ! else + ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) + ! ! node => list_next(previous_node) + ! end if + cont_db(3) = cont_db(3) + 6 + cont_int(3) = cont_int(3) + 4 + end if + if (cell(2) <= domx(1) .and. cell(2) > 1) then + ! print*, "L 554", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + ! print*, "domx", domx + LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1), x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "id",id,"transferindo para o oeste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! if (cell(2) < domx(1)) then + ! ! print*, 'affs' + ! ! read(*,*) + ! ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_W + + ! ! call list_remove(previous_node) + ! ! node => list_next(previous_node) + ! ! else + ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) + ! ! node => list_next(previous_node) + ! end if + cont_db(4) = cont_db(4) + 6 + cont_int(4) = cont_int(4) + 4 + end if + if (cell(1) >= domy(2) .and. cell(1) < mesh(2)+2) then + ! print*, "L 567", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + ! print*, "id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! if (cell(1) > domy(1)) then + ! ! read(*,*) + ! ! print*, 'affs' + ! ! call list_remove(previous_node) + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_N(cont_int(1)+1:cont_int(1)+4), "cont", cont_int(1) + ! ! node => list_next(previous_node) + ! ! else + ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) + ! ! node => list_next(previous_node) + ! end if + cont_db(1) = cont_db(1) + 6 + cont_int(1) = cont_int(1) + 4 + end if + if (cell(1) <= domy(1) .and. cell(1) > 1) then + ! print*, "L 580", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + LT%lstrdb_S(cont_db(2)+1:cont_db(2)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_S(cont_int(2)+1:cont_int(2)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L 583" + ! print*, "id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! if (cell(1) < domy(1)) then + ! ! read(*,*) + ! ! print*, 'affs' + ! ! call list_remove(previous_node) + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_S(cont_int(2)+1:cont_int(2)+4), "cont", cont_int(1) + ! ! node => list_next(previous_node) + ! ! else + ! ! !!!! ! print*, "L 590" + ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) + ! ! node => list_next(previous_node) + ! ! !!!! ! print*, "L 624" + ! end if + cont_db(2) = cont_db(2) + 6 + cont_int(2) = cont_int(2) + 4 + end if + + ! FIM DA VERIFICAÇÃO SE MUDOU DE DOMÍNIO + + ! Antes aqui tinha um else para o caso de não mudar de processo. + ! Removi todos os call list_remove e deixei que a partícula fosse para + ! a área emprestada que será limpa na próxima execução de call comp_x. Isto pois a célula precisa + ! sentir a presença da partícula que ela perdeu mas ficou ao lado + ! print*, "mudou2" + ! if (id == 0) read(*,*) + + call list_change(previous_node,malha(cell(1),cell(2))%list) + node => list_next(previous_node) + + else !Não mudou de celula + ! Mas tem que avisar quais ainda estão na + ! vizinhança + + ! print*, "NAO mudou cell",cell, "domx", domx, "domy", domy + if (cell(2) == domx(2)) then + !!! ! print*, "L 538", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + ! 6 elementos + LT%lstrdb_E(cont_db(3)+1:cont_db(3)+6) = [x(1),x(2), & + ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + ! 4 elementos + LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(3) = cont_db(3) + 6 + cont_int(3) = cont_int(3) + 4 + + else if (cell(2) == domx(1)) then + !!! ! print*, "L 554", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o oeste", LT%lstrint_W + cont_db(4) = cont_db(4) + 6 + cont_int(4) = cont_int(4) + 4 + end if + + if (cell(1) == domy(1)) then + !!! ! print*, "L 567", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_S(cont_db(2)+1:cont_db(2)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_S(cont_int(2)+1:cont_int(2)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(2) = cont_db(2) + 6 + cont_int(2) = cont_int(2) + 4 + + elseif (cell(1) == domy(2)) then + !!! ! print*, "L 580", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(1) = cont_db(1) + 6 + cont_int(1) = cont_int(1) + 4 + end if + + previous_node => node + node => list_next(node) + ! print*, "L740 FORÇA", ptr%p%F, id + + end if + if (associated(node)) then + ptr = transfer(list_get(node), ptr) + ! aqui vemos se já passou o tempo que as posições ficaram travadas. + if (propriedade(ptr%p%grupo)%x_lockdelay > t) then + ptr%p%flag = .false. + end if + end if + end do + end do + end do + + ! transferir os termos que estão no encontro de 4 processos + ! print*, "L 794", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + if (np > 1 .and. sum(ids(5:8)) > -4) then + + do i = 5,8 + !identifica qual o processo na diagonal + if (ids(i) > -1) then + destD = i + end if + end do + ! ! print*, "L 787", ids, "D", destD + ! Aqui está com suporte a 4 processadores por enquanto + ! Primeiro transfere para as celulas emprestadas + + if (domy(1) /= 1) i = domy(1) + if (domy(2) /= mesh(2)+2) i = domy(2) + if (domx(1) /= 1) j = domx(1) + if (domx(2) /= mesh(1) + 2) j = domx(2) + + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + + if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "transferencia diagonal", i,j, "ID", id + do while(associated(node)) + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + ! Segundo transfere para as celulas do domínio caso a partícula tenha mudado de processo + if (domy(1) /= 1) i = domy(1) -1 + if (domy(2) /= mesh(2)+2) i = domy(2) +1 + if (domx(1) /= 1) j = domx(1) -1 + if (domx(2) /= mesh(1) + 2) j = domx(2) +1 + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "mudança diagonal", i,j, "ID", id + do while(associated(node)) + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"mudando diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + ! Fim da trasnfefência para diagonal + end if + + ! print*, "L 854", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, "L SOMA", (ids(1) + ids(2) + ids(3) +ids(4)), "id", id + ! print*, "enviando" + if (np > 1) then !se for paralelo + ! call MPI_SEND(BUF, COUNT, DATATYPE, DEST, TAG, COMM, IERROR) + tag = 1 + if (ids(1) >= 0) call MPI_SEND(LT%lstrdb_N, cont_db(1), MPI_DOUBLE_PRECISION, ids(1), tag,MPI_COMM_WORLD, ierr) + if (ids(2) >= 0) call MPI_SEND(LT%lstrdb_S, cont_db(2), MPI_DOUBLE_PRECISION, ids(2), tag,MPI_COMM_WORLD, ierr) + if (ids(3) >= 0) call MPI_SEND(LT%lstrdb_E, cont_db(3), MPI_DOUBLE_PRECISION, ids(3), tag,MPI_COMM_WORLD, ierr) + if (ids(4) >= 0) call MPI_SEND(LT%lstrdb_W, cont_db(4), MPI_DOUBLE_PRECISION, ids(4), tag,MPI_COMM_WORLD, ierr) + if (sum(ids(5:8)) > -4) then + call MPI_SEND(LT%lstrdb_D, cont_db(destD), MPI_DOUBLE_PRECISION, ids(destD), tag,MPI_COMM_WORLD, ierr) + end if + tag = 2 + if (ids(1) >= 0) call MPI_SEND(LT%lstrint_N, cont_int(1), MPI_integer, ids(1), tag,MPI_COMM_WORLD, ierr) + if (ids(2) >= 0) call MPI_SEND(LT%lstrint_S, cont_int(2), MPI_integer, ids(2), tag,MPI_COMM_WORLD, ierr) + if (ids(3) >= 0) call MPI_SEND(LT%lstrint_E, cont_int(3), MPI_integer, ids(3), tag,MPI_COMM_WORLD, ierr) + if (ids(4) >= 0) call MPI_SEND(LT%lstrint_W, cont_int(4), MPI_integer, ids(4), tag,MPI_COMM_WORLD, ierr) + if (sum(ids(5:8)) > -4) then + call MPI_SEND(LT%lstrint_D, cont_int(destD), MPI_integer, ids(destD), tag,MPI_COMM_WORLD, ierr) + end if + + ! print*, "L 935 TUDO ENVIADO", id + tag = 1 + if (ids(1) >= 0) then + call MPI_probe(ids(1),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(1) , ierr ) + call MPI_RECV (LT%lstrdb_N, cont_db(1), MPI_DOUBLE_PRECISION,ids(1), tag, MPI_COMM_WORLD, status, ierr) + end if + + if (ids(2) >= 0) then + call MPI_probe(ids(2),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(2) , ierr ) + call MPI_RECV (LT%lstrdb_S, cont_db(2), MPI_DOUBLE_PRECISION,ids(2), tag, MPI_COMM_WORLD, status, ierr) + + end if + if (ids(3) >= 0) then + call MPI_probe(ids(3),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(3) , ierr ) + call MPI_RECV (LT%lstrdb_E, cont_db(3), MPI_DOUBLE_PRECISION,ids(3), tag, MPI_COMM_WORLD, status, ierr) + end if + + if (ids(4) >= 0) then + call MPI_probe(ids(4),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(4) , ierr ) + call MPI_RECV (LT%lstrdb_W, cont_db(4), MPI_DOUBLE_PRECISION,ids(4), tag, MPI_COMM_WORLD, status, ierr) + end if + + tag = 2 + if (ids(1) >= 0) then + call MPI_probe(ids(1),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(1) , ierr ) + call MPI_RECV (LT%lstrint_N, cont_int(1), MPI_integer,ids(1), tag, MPI_COMM_WORLD, status, ierr) + end if + if (ids(2) >= 0) then + call MPI_probe(ids(2),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(2) , ierr ) + call MPI_RECV (LT%lstrint_S, cont_int(2), MPI_integer,ids(2), tag, MPI_COMM_WORLD, status, ierr) + end if + if (ids(3) >= 0) then + call MPI_probe(ids(3),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(3) , ierr ) + call MPI_RECV (LT%lstrint_E, cont_int(3), MPI_integer,ids(3), tag, MPI_COMM_WORLD, status, ierr) + end if + if (ids(4) >= 0) then + call MPI_probe(ids(4),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(4) , ierr ) + call MPI_RECV (LT%lstrint_W, cont_int(4), MPI_integer,ids(4), tag, MPI_COMM_WORLD, status, ierr) + end if + if (sum(ids(5:8)) > -4) then + tag = 1 + ! RECEBE DA DIAGONAL (suporte para 4 processos) + ! call MPI_probe(,tag, MPI_COMM_WORLD, status, ierr) + call MPI_probe(ids(destD),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(destD) , ierr ) + call MPI_RECV (LT%lstrdb_D, cont_db(destD), MPI_DOUBLE_PRECISION,ids(destD), tag, MPI_COMM_WORLD, status, ierr) + tag = 2 + ! Recebe INT da DIAGONAL. Aqui está com suporte para 4 processadores + call MPI_probe(ids(destD),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_int(destD) , ierr ) + call MPI_RECV (LT%lstrint_D, cont_int(destD), MPI_DOUBLE_PRECISION,ids(destD), tag, MPI_COMM_WORLD, status, ierr) + + ! DIAGONAL, tem que alterar para o caso com >4 processos + j = destD + do i = 0, cont_db(j)/6 + if (cont_db(j) > 0) then + ! print*, "LA 1108" + ! print*, "D al", id, LT%lstrdb_D(i*6+1:i*6+2) + allocate(ptr%p) + ptr%p%x = LT%lstrdb_D(i*6+1:i*6+2) + ! print*, ptr%p%x + ptr%p%v = LT%lstrdb_D(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_D(i*4+4) + ptr%p%F = LT%lstrdb_D(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_D(i*4+3) !identidade da partícula importante para imprimir + call list_insert(malha(LT%lstrint_D(i*4+1), & + LT%lstrint_D(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "D allocado em", LT%lstrint_D(i*4+1), LT%lstrint_D(i*4+2), "id =", id, "posição",ptr%p%x + ! print*, "W allocado F = ",ptr%p%F + cont_db(j) = cont_db(j) - 6 + end if + end do + end if + + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, 'ID APOS BARREIRA', id + ! print*, "L 1033" + do j = 1,4 + do i = 0, cont_db(j)/6 + !! ! print*, "L 840" , i, j, "id", id, "cont_db/6", cont_db(j)/6 + if (j == 1 .and. cont_db(1) > 0 ) then + ! print*, "LA 1133" + allocate(ptr%p) + + ptr%p%x = LT%lstrdb_N(i*6+1:i*6+2) + ptr%p%v = LT%lstrdb_N(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_N(i*4+4) + ptr%p%F = LT%lstrdb_N(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_N(i*4+3) !identidade da partícula importante para imprimir + ! !print*, 'LLLL', LT%lstrint_N(i*4+1), LT%lstrint_N(i*4+2) + call list_insert(malha(LT%lstrint_N(i*4+1), & + LT%lstrint_N(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "N allocado em", LT%lstrint_N(i*4+1), LT%lstrint_N(i*4+2), "id =", id + ! print*, "N allocado F = ",ptr%p%F + cont_db(j) = cont_db(j) - 6 + + end if + if (j == 2 .and. cont_db(2) > 0) then + ! print*, "LA 1151" + allocate(ptr%p) + !! ! print*, "L 859 <<<" + ptr%p%x = LT%lstrdb_S(i*6+1:i*6+2) + ptr%p%v = LT%lstrdb_S(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_S(i*4+4) + ptr%p%F = LT%lstrdb_S(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_S(i*4+3) !identidade da partícula importante para imprimir + !! ! print*, "L 865 <<<", LT%lstrint_S + call list_insert(malha(LT%lstrint_S(i*4+1), & + LT%lstrint_S(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "S allocado em", LT%lstrint_S(i*4+1), LT%lstrint_S(i*4+2), "id =", id,i + ! print*, "S allocado F = ",ptr%p%F + !! ! print*, "L 869 <<<" + cont_db(j) = cont_db(j) - 6 + end if + if (j == 3 .and. cont_db(3) > 0) then + ! print*, "LA 1170" + allocate(ptr%p) + ptr%p%x = LT%lstrdb_E(i*6+1:i*6+2) + ptr%p%v = LT%lstrdb_E(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_E(i*4+4) + ptr%p%F = LT%lstrdb_E(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_E(i*4+3) !identidade da partícula importante para imprimir + call list_insert(malha(LT%lstrint_E(i*4+1), & + LT%lstrint_E(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "E allocado em", LT%lstrint_E(i*4+1), LT%lstrint_E(i*4+2), "x=",ptr%p%x + ! print*, "E allocado", LT%lstrint_E + cont_db(j) = cont_db(j) - 6 + end if + if (j == 4 .and. cont_db(4) > 0 ) then + ! print*, "LA 1185" + allocate(ptr%p) + ptr%p%x = LT%lstrdb_W(i*6+1:i*6+2) + ! print*, ptr%p%x + ptr%p%v = LT%lstrdb_W(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_W(i*4+4) + ptr%p%F = LT%lstrdb_W(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_W(i*4+3) !identidade da partícula importante para imprimir + call list_insert(malha(LT%lstrint_W(i*4+1), & + LT%lstrint_W(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "W allocado em", LT%lstrint_W(i*4+1), LT%lstrint_W(i*4+2), "id =", id + ! print*, "W allocado F = ",ptr%p%F + cont_db(j) = cont_db(j) - 6 + end if + end do + end do + ! print*, "L 1106" + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end if + + ! print*, "L 1115 fim comp_x", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end subroutine comp_x + + ! atualiza velocidades + subroutine comp_v(malha,mesh,dt,t,propriedade,domx,domy) + use linkedlist + use mod1 + use data + use mod0 + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer :: i,j + integer, intent(in) :: mesh(2),domx(2),domy(2) + type(list_t), pointer :: node + real(dp), intent(in) :: dt, t + type(data_ptr) :: ptr + + do i = domy(1),domy(2) + do j = domx(1),domx(2) + 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 + m = propriedade(ptr%p%grupo)%m + ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) + ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) + end if + ptr%p%F = [0,0] + node => list_next(node) + end do + end do + end do + end subroutine comp_v + + subroutine walls(icell,jcell,mesh,malha,domx,domy,north,south,east,west,id) + use linkedlist + use mod1 + use data + use mod0 + + type(container), allocatable,dimension(:,:),intent(in) :: malha + real(dp),intent(in) :: icell(:), jcell(:) + character(1), intent(in) :: north, south, east, west + type(list_t), pointer :: node, previous_node + integer :: i,j,cell(2) + integer, intent(in) :: mesh(:) + type(data_ptr) :: ptr,ptrn + integer, intent(in) :: domx(2), domy(2), id + + ! ! print*, "L 1004",id + if (domy(2) == mesh(2)+2) then + !! ! print*, "L 980", id + if (north == 'e') then !elastic + i = mesh(2) +2 + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + + do while (associated(node)) + ! print*, "cell",i,j, "north" + ptr = transfer(list_get(node), ptr) + ptr%p%x(2) = 2*icell(mesh(2)+1) - ptr%p%x(2) + ptr%p%v(2) = -ptr%p%v(2) + call list_change(previous_node,malha(i-1,j)%list) + node => list_next(previous_node) + ! print*, 'part=',ptr%p%n, "cell",i,j + end do + end do + else if (north == 'o') then + i = mesh(2) +2 + do j = domx(1), domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + do while (associated(node)) + ! print*, "cell",i,j + ptr = transfer(list_get(node), ptr) + ptr%p%x(2) = 0.0/0.0 + ptr%p%v(2) = 0.0/0.0 + call list_remove(previous_node) + node => list_next(previous_node) + end do + end do + else !periodic + print*, oi + end if + end if + + !! print*, "Linha 1018",id + if (domy(1) == 1) then + !! ! print*, "L 1019", id + if (south == 'e') then + i = 1 + do j = domx(1), domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + do while (associated(node)) + ! print*, "cell",i,j, "south", id + ptr = transfer(list_get(node), ptr) + ! print*, "cell",i,j + ! print*, 'part=',ptr%p%n + + ptr%p%x(2) = -ptr%p%x(2) + ptr%p%v(2) = -ptr%p%v(2) + call list_change(previous_node,malha(2,j)%list) + + node => list_next(previous_node) + end do + end do + else if (south == 'o') then + i = 1 + do j = domx(1), domx(2) + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + ptr%p%x(2) = 0.0/0.0 + ptr%p%v(2) = 0.0/0.0 + call list_remove(previous_node) + node => list_next(previous_node) + end do + end do + else !periodic + print*, oi + end if + end if + !! ! print*, "L 1057", id + if (domx(1) == 1) then + if (west == 'e') then + j = 1 + do i = domy(1), domy(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + ! print*, "n", ptr%p%x + ptr%p%x(1) = -ptr%p%x(1) + ptr%p%v(1) = -ptr%p%v(1) + call list_change(previous_node,malha(i,2)%list) + node => list_next(previous_node) + end do + end do + else if (west == 'o') then + j = 1 + do i = domy(1), domy(2) + previous_node => malha(i,j)%list + + node => list_next(malha(i,j)%list) + ! ! !print*, 'Linha 236', associated(previous_node) + do while (associated(node)) + ! print*,'N'!,x(ptr,2) + !read(*,*) + ptr = transfer(list_get(node), ptr) + ptr%p%x(2) = 0.0/0.0 + ptr%p%v(2) = 0.0/0.0 + call list_remove(previous_node) + ! previous_node => node + node => list_next(previous_node) + end do + end do + else !periodic + print*, oi + end if + end if + if (domx(2) == mesh(1) +2) then + if (east == 'e') then + j = mesh(1)+2 + do i = domy(1), domy(2) + ! print*, "A" + previous_node => malha(i,j)%list + ! print*, "B" + node => list_next(malha(i,j)%list) + ! print*, "C" + do while (associated(node)) + ! print*, "cell",i,j, "east" + ptr = transfer(list_get(node), ptr) + ! print*, "E" + ptr%p%x(1) = 2*jcell(mesh(1)+1) - ptr%p%x(1) + ptr%p%v(1) = -ptr%p%v(1) + ! print*, "F" + call list_change(previous_node,malha(i,j-1)%list) + ! print*, "G" + node => list_next(previous_node) + end do + end do + else if (east == 'o') then + j = mesh(1)+2 + ! print*,'j = ',j + do i = domy(1), domy(2) + previous_node => malha(i,j)%list + + node => list_next(malha(i,j)%list) + !print*, 'Linha 236', associated(previous_node) + do while (associated(node)) + ! print*,'N'!,x(ptr,2) + !read(*,*) + ptr = transfer(list_get(node), ptr) + ptr%p%x(1) = 0.0/0.0 + ptr%p%v(1) = 0.0/0.0 + call list_remove(previous_node) + ! previous_node => node + node => list_next(previous_node) + end do + end do + else !periodic + print*, oi + end if + end if + ! print*,"linha 1146 ok" + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! !read(*,*) + end subroutine walls +! + subroutine distribuicao_inicial(malha,N,mesh) + use linkedlist + use mod1 + use data + use mod0 + + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer :: i,j,aux1 + integer, intent(in) :: N,mesh(2) + type(list_t), pointer :: node + type(data_ptr) :: ptr + + + do i = 1,mesh(2)+2 + do j = 1,mesh(1)+2 + node => list_next(malha(i,j)%list) + aux1 = 0 + do while (associated(node)) + aux1 = aux1 +1 + ptr = transfer(list_get(node), ptr) + ! celula(ptr,:) = [j,i] + print *, 'posição', i, ',', j + print *, ptr%p%x + node => list_next(node) + if (aux1 > N) then + print*, 'deu ruim' + exit + end if + end do + end do + end do + end subroutine distribuicao_inicial + + subroutine clean_mesh(malha, mesh, domx, domy,id, tudo) + ! This subroutine cleans the unaccessible cells of each process bounded by domx, domy + use linkedlist + use mod1 + use data + use mod0 + + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer :: i,j, k = 0 + integer, intent(in) :: mesh(2),domx(2),domy(2) + type(list_t), pointer :: node, previous_node + type(data_ptr) :: ptr + logical, intent(in) :: tudo + ! print*, "AAA", id + if (tudo) k = 1 + do i = 2-k,mesh(1)+1+k + do j = 2-k, mesh(2)+1+k + if (((i < domy(1)-1 .or. i > domy(2)+1) .and. (j < domx(1)-1 .or. j > domx(2)+1)) .or. tudo) then + node => list_next(malha(i,j)%list) + ! print*, "i,j,id", i,j, id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + node => list_next(node) + end do + node => list_next(malha(i,j)%list) + if (associated(node)) then + call list_free(malha(i,j)%list) + call list_init(malha(i,j)%list) + ! if (.not. tudo) then + ! print*, "BBB" + ! else + ! print*, "CCC" + ! end if + end if + end if + end do + end do + + if (tudo) then + do i = 2-k,mesh(1)+1+k + do j = 2-k, mesh(2)+1+k + deallocate(malha(i,j)%list) + end do + end do + ! dellocate( ) + end if + + end subroutine clean_mesh + +end module mod2 + + +! Este módulo escreve csv +module saida + contains + subroutine vec2csv(v,n,d,prop,step,t,nimpre,start) + use mod1 + !N é o número de partículas e d é a dimensão da propriedade + ! cada linha é uma partícula + implicit none + integer, intent(in) :: n,d + real(dp), intent(in) :: v(n,d), t, start + integer :: i,step, nimpre + character(*) :: prop + character(4) :: extensao = '.csv', passo + real(dp) :: time + real(dp), save :: timep, etc, dtimepp + + write(passo,'(i0)') step + open(10,file='temp/'//prop//extensao//'.'//trim(passo),status="replace") + if (d == 1) then + do i = 1,n + write(10,*) v(i,1) + end do + else if (d == 2) then + do i = 1,n + write(10,*) v(i,1),',',v(i,2) + end do + else + do i = 1,n + write(10,*) v(i,1),',',v(i,2),',',v(i,3) + end do + end if + close(10) + + if (prop == "position") then + if (step == 0) timep = 0 + call cpu_time(time) + if (step == 1) then + etc = ((time - start)/step+ (time-timep))*0.5*nimpre - (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 + etc = ((time - start)*2/step + ((time-timep)*4 + dtimepp*4))*(nimpre-step)/10 + print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", f10.3, "s" )',prop//extensao//'.'//trim(passo),t,etc + dtimepp = (time-timep) + else + print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", "unknown" )',prop//extensao//'.'//trim(passo),t + end if + + timep = time + else + if (step > 0) then + print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", f10.3, "s" )',prop//extensao//'.'//trim(passo),t,etc + else + print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", "unknown" )',prop//extensao//'.'//trim(passo),t + end if + timep = time + end if + ! print*, 'Salvo arquivo ',prop//extensao//'.'//trim(passo), "t =", t, "ETC: ", etc + + end subroutine vec2csv + + subroutine linked2vec(malha,domx,domy,nxv,aux1) + + use linkedlist + use mod1 + use data + use mod0 + + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer :: i,j,aux1 + integer, intent(in) :: domx(2), domy(2) + type(list_t), pointer :: node + real(dp), intent(out) :: nxv(:) + type(data_ptr) :: ptr + + aux1 = 1 +! real(dp),intent(inout) :: celula(:,:) + + do i = domy(1), domy(2) + do j = domx(1), domx(2) + ! print *, 'posição', i, ',', j + node => list_next(malha(i,j)%list) + ! print*, i,j + ! print*, 'is associated?', associated(node) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + nxv(aux1:aux1+4) = [ real(ptr%p%n, kind(0.d0)), ptr%p%x(1),ptr%p%x(2), & + ptr%p%v(1), ptr%p%v(2)] + aux1 = aux1 + 5 + node => list_next(node) + end do + end do + end do + aux1 = aux1 -1 + + end subroutine linked2vec + + subroutine nancheck(malha,mesh,t) + use linkedlist + use mod1 + use data + use mod0 + + real(dp), intent(in) :: t + integer :: i + integer, intent(in) :: mesh(:) + type(container), allocatable,dimension(:,:),intent(in) :: malha + type(data_ptr) :: ptr + type(list_t), pointer :: node + do i = 1,mesh(2)+2 + do j = 1,mesh(1)+2 + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + !calcula a energia cinética atual + node => list_next(node) + if (isnan(ptr%p%x(1))) then + print*, 'NaN em part',ptr%p%n,'x, t=',t + end if + if (isnan(ptr%p%x(2))) then + print*, 'NaN em part',ptr%p%n,'y, t=',t + end if + end do + end do + end do + + + + do i = 1,N + + end do + end subroutine nancheck + +end module saida + +program main + use mod1 + use linkedlist + use data + use mod0 + use matprint + use saida + use m_config + use mod2 + use randnormal + use mpi + + 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) + integer :: subx, suby, NMPT, j2 + integer, target :: k + real(dp), dimension(:,:), allocatable :: v, x, celula !força n e n+1, velocidades e posições + real(dp), dimension(:), allocatable :: icell,jcell, nxv, nxv_send !dimensões das celulas e vetor de resultado pra imprimir + real(dp) :: t=0,t_fim,dt,printstep, sigma, epsil, rcut,aux2,start = 0,finish,Td,kb = 1.38064852E-23,fric_term,vd(2) + real(dp) :: GField(2), temp_Td(3), dimX, dimY + integer,allocatable :: interv(:), interv_Td(:), grupo(:), rcounts(:), displs(:) + type(container), allocatable,dimension(:,:) :: malha + type(prop_grupo),allocatable,dimension(:) :: propriedade + type(list_t), pointer :: node + ! Variáveis do arquivo de configuração + type(CFG_t) :: my_cfg + character(5) :: particle + character(4) :: wall + character(11) :: nome,arquivo !nome de partícula deve ter até 11 caracteres + type(string) :: part_nomes(10) ! vetor de strings + character(1) :: optio + type(data_ptr) :: ptr !integer,pointer :: ptr,ptrn + logical :: laux, termostato + character(len=32) :: arg + character :: arg1(1) + type(lstr) :: LT + integer ( kind = 4 ) status(MPI_STATUS_SIZE) + integer ( kind = 4 ) ierr + integer ( kind = 4 ) np + integer ( kind = 4 ) tag + integer ( kind = 4 ) id, ids(8) + + ! depois ver se não da pra fazer tudo serialmente + call MPI_Init ( ierr ) + call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr ) + call MPI_Comm_size ( MPI_COMM_WORLD, np, ierr ) + + ! argumento para checar NAN no nancheck + call get_command_argument(1, arg) + arg1 = trim(arg) + +! CONDIÇÕES INICIAIS/CORTORNO E PROPRIEDADES + ! Propriedades que cobrem todo um grupo de partículas + ! serão armazenadas na estrutura propriedade(grupo) +! Inicia o método de configuração + call CFG_add(my_cfg, "global%N",1,& + "number of particles") + call CFG_add(my_cfg, "global%Ntype",1,& + "number of types") + call CFG_add(my_cfg, "global%dt",0.0_dp,& + "time step") + call CFG_add(my_cfg, "global%t_fim",0.0_dp,& + "simulation time") + call CFG_add(my_cfg, "global%nimpre",1,& + "number of result files") + call CFG_add(my_cfg, "global%dimX",1.0_dp,& + "y-dimension of the region of calculus") + call CFG_add(my_cfg, "global%dimY",1.0_dp,& + "y-dimension of the region of calculus") + call CFG_add(my_cfg, "global%mesh",(/0, 0/),& + "mesh size") !celulas + call CFG_add(my_cfg, "global%sigma",1.0_dp,& + "Leonard-Jones sigma factor") + call CFG_add(my_cfg, "global%epsilon",1.0_dp,& + "Leonard-Jones epsilon factor") + call CFG_add(my_cfg, "global%rcut",1.0_dp,& + "Cut radius") + call CFG_add(my_cfg,"global%wall",'eeee', & + "wall's periodic vs elastic") + call CFG_add(my_cfg,"global%Td",-1.0_dp, & + "Thermostat") + call CFG_add(my_cfg,"global%temp_Td",(/0.0_dp, 0.0_dp, 0.0_dp/), & + "Thermostat iteractions time") + call CFG_add(my_cfg,"global%vd",(/0.0_dp, 0.0_dp/), & + "Mean vd") + call CFG_add(my_cfg,"global%fric_term",1.0_dp, & + "Friction term") + call CFG_add(my_cfg,"global%NMPT",1, & + "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_read_file(my_cfg, "settings.ini") + call CFG_get(my_cfg,"global%N",N) + call CFG_get(my_cfg,"global%Ntype",Ntype) + call CFG_get(my_cfg,"global%t_fim",t_fim) + call CFG_get(my_cfg,"global%dt",dt) + call CFG_get(my_cfg,"global%nimpre",nimpre) + call CFG_get(my_cfg,"global%dimX", dimX) + call CFG_get(my_cfg,"global%dimY", dimY) + call CFG_get(my_cfg,"global%mesh",mesh) + call CFG_get(my_cfg,"global%rcut",rcut) + call CFG_get(my_cfg,"global%wall",wall) + call CFG_get(my_cfg,"global%Td",Td) + call CFG_get(my_cfg,"global%temp_Td",temp_Td) + call CFG_get(my_cfg,"global%vd",vd) + call CFG_get(my_cfg,"global%fric_term",fric_term) + call CFG_get(my_cfg,"global%NMPT",NMPT) + call CFG_get(my_cfg,"global%GField",GField) + + + if (NMPT < 1) NMPT = N + printstep = ((temp_Td(2)-temp_Td(1))/dt)/temp_Td(3) + if (printstep < 1) printstep = 1 + !aloca as Variáveis com base no número de partículas + allocate(v(N,2),x(N,2),grupo(N),interv(nimpre),interv_Td(nint(printstep)), & + jcell(mesh(2)+1),icell(mesh(1)+1),celula(N,2),propriedade(Ntype)) + allocate(LT%lstrdb_N(NMPT*6),LT%lstrdb_S(NMPT*6),LT%lstrdb_E(NMPT*6),LT%lstrdb_W(NMPT*6),LT%lstrdb_D(NMPT*6), & + LT%lstrint_N(NMPT*4),LT%lstrint_S(NMPT*4),LT%lstrint_E(NMPT*4),LT%lstrint_W(NMPT*4), LT%lstrint_D(NMPT*4)) + ! allocate(nxv(N*5),nxv_send(N*5)) + + aux2 = dimX/mesh(1) + jcell = (/((i*aux2),i = 0,mesh(1))/) ! direção x + + aux2 = dimY/mesh(2) + icell = (/((i*aux2),i = 0,mesh(2))/) ! direção Y + + + do i = 0,(Ntype-1) + if (id == 0) write(*,'(a,i0)') 'par_',i + write(particle,'(a,i0)') 'par_',i + call CFG_add(my_cfg, particle//"%m",1.1_dp,& + "mass of "//particle) + call CFG_add(my_cfg, particle//"%v",(/0.0_dp, 0.0_dp/), & + "velocity of "//particle) + call CFG_add(my_cfg, particle//"%x","arquivo.csv", & + "position of "//particle) + call CFG_add(my_cfg, particle//"%nome",'abcde', & + "name of "//particle) + call CFG_add(my_cfg, particle//"%quantidade",1, & + "quantity of "//particle) + call CFG_add(my_cfg, particle//"%epsilon",1.1_dp, & + "sigma "//particle) + call CFG_add(my_cfg, particle//"%sigma",1.1_dp, & + "epsilon"//particle) + call CFG_add(my_cfg, particle//"%x_lockdelay",1.1_dp, & + "change in position delay "//particle) + end do + + + 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 + call CFG_get(my_cfg, particle//"%x", arquivo) + call CFG_get(my_cfg, particle//"%m", propriedade(i+1)%m) + call CFG_get(my_cfg, particle//"%epsilon", propriedade(i+1)%epsilon) + call CFG_get(my_cfg, particle//"%sigma", propriedade(i+1)%sigma) + call CFG_get(my_cfg, particle//"%x_lockdelay", propriedade(i+1)%x_lockdelay) + ! le o arquivo com posições + + open(20,file=arquivo,status='old') + do j = 1,quant + call CFG_get(my_cfg, particle//"%v", v(cont,:)) + grupo(cont) = i+1 + read(20,*) x(cont,:) !,x(cont,2),x(cont,3) + cont = cont+1 + end do + end do + + ! mostra informações lidas + if (id == 0) then + call CFG_write(my_cfg, "settings.txt") + print*, "Nomes das partículas" + do i = 1,Ntype + print*, part_nomes(i)%str + end do + end if + + !--------------------------------------------------! + ! CRIA LISTAS LIGADAS + + !cada elemento da matriz representa uma celula da malha + allocate(malha(mesh(2)+2,mesh(1)+2)) ! o primeiro e últimos índices correspondem a + ! celulas fantasma + !inicialização + k = 0 + do i = 1,mesh(2)+2 !linha + do j = 1,mesh(1)+2 + ! allocate(ptr%p) + ! ptr%p%x = [0,0] + ! ptr%p%v = [0,0] + ! ptr%p%grupo = 0 + ! ptr%p%F = [0,0] + ! ptr%p%n = k + ! ptr%p%flag = .true. ! flag auxiliar para uso + call list_init(malha(i,j)%list) !, DATA=transfer(ptr, list_data)) + end do + end do + ! Vetor que indica quando imprimir + printstep = nint(t_fim/dt)/nimpre + interv = (/(printstep*i,i=0,nimpre) /) + ! Vetor que indica quando fazer velocity scaling + printstep = nint(temp_Td(3)) + interv_Td = (/(printstep*i,i=0,nint( ((temp_Td(2)-temp_Td(1))/dt)/temp_Td(3) )) /) + interv_Td = interv_Td + nint(temp_Td(1)/dt) + i = 0 + + if (Td == 0) then + Td = (2/(3*N*kb))*sum(0.5*propriedade(ptr%p%grupo)%m*(v(:,1)**2+v(:,2)**2)) + else if (Td < 0) then + interv_Td = -1 + end if +! print*, "2073" + !Aloca as partículas nas celulas + do k = 1,cont-1 + + laux = .true. + + j = 2 + do while (laux) + if (x(k,1) > jcell(j)) then + j = j+1 + else + laux = .false. + end if + end do + + laux = .true. + i = 2 + do while (laux) + if (x(k,2) > icell(i)) then + i = i+1 + else + laux = .false. + end if + end do + !i e j contém agora a posição da partícula na matriz malha + + allocate(ptr%p) + ptr%p%x = x(k,:) + ptr%p%v = v(k,:) + ptr%p%grupo = grupo(k) + ptr%p%F = [0,0] !passo 1 do SV + ptr%p%n = k !identidade da partícula importante para imprimir + call list_insert(malha(i,j)%list, data=transfer(ptr, list_data)) + end do + + deallocate(grupo) + + ! adiciona às particulas velocidade inicial de acordo com distribuição maxwell boltzmann + if (vd(1) /= 0 .and. vd(2) /= 0 .and. vd(1) /= 0) call MaxwellBoltzmann(malha,mesh,sqrt(vd)) + + !-------------------------------------------------! + ! ITERAÇÕES + + ! imprime condições iniciais + + + if (id == 0) then + + j = 0 + call system('mkdir temp') !pasta temporária para armazenar os resultados + ! call linked2vec(malha,domx,domy,nxv,aux1) + call vec2csv(x,N,2,'position',j,t,nimpre,start) + call vec2csv(v,N,2,'velocity',j,t,nimpre,start) + ! call vec2csv(celula,N,2,'cell',j) + j = j + 1 + + ! print*,'V_tot0 = ',comp_pot(mesh,malha,propriedade,rcut) + end if + + !! AQUI COMEÇA O PARALELISMO + + ! só vai funcionar com serial e np par + ! Descobir a melhor forma de dividir a malha + + if (dimX > 3*dimY .and. np > 1) then + subx = 4 + suby = 1 + else if (dimY > 3*dimX .and. np > 1) then + subx = 1 + suby = 4 + else if (np > 1) then + subx = 2 + suby = 2 + else + subx = 1 + suby = 1 + end if + + np = subx*suby + ! print*, "NP", np + ids = [-1,-1,-1,-1,-1,-1,-1,-1] + + if (np > 2) then + L_O: do j = 0, (suby-1) + do i = 0, (subx-1) + if (id == (i+j*suby)) then + ! ids = [N,S,E,W, NE, NW, SE, SW] + ! Depois tentar para caso geral + ! ids = [(j-1)*subx +i, (j+1)*subx +i, & + ! j*subx + i +1, -(i>0)*(j*subx + i)-1] + ! if (ids(1) < 0) ids(1) = -1 + ! if (ids(4) < 0) ids(4) = -1 + ! if (j == suby-1) ids(2) = -1 + ! if (i == subx -1) ids(3) = -1 + if (subx == 2) then + if (id == 0) then + ids = [2, -1, 1,-1, 3, -1, -1, -1] + else if (id == 1) then + ids = [3, -1, -1, 0, -1, 2, -1, -1] + else if (id == 2) then + ids = [-1, 0, 3, -1, -1, -1, 1, -1] + else + ids = [-1, 1, -1, 2, -1, -1, -1, 0] + end if + else if (subx == 4) then + if (id == 0) then + ids = [-1, -1, 1,-1, -1, -1, -1, -1] + else if (id == 1) then + ids = [-1, -1, 2, 0, -1, -1, -1, -1] + else if (id == 2) then + ids = [-1, -1, 3, 1, -1, -1, -1, -1] + else + ids = [-1, -1, -1, 2, -1, -1, -1, -1] + end if + else if (suby == 4) then + if (id == 0) then + ids = [1, -1, -1, -1, -1, -1, -1, -1] + else if (id == 1) then + ids = [2, 0, -1, -1, -1, -1, -1, -1] + else if (id == 2) then + ids = [3, 1, -1, -1, -1, -1, -1, -1] + else + ids = [-1, 2, -1, -1, -1, -1, -1, -1] + end if + end if + domx = [i*int((mesh(1)+2)/subx)+1, (i+1)*int((mesh(1)+2)/(subx))] + domy = [j*int((mesh(2)+2)/suby)+1, (j+1)*int((mesh(2)+2)/(suby))] + exit L_O + end if + end do + end do L_O + else + domx = [1,mesh(1)+2] + domy = [1,mesh(2)+2] + end if + + print '("id = ",I2, "; domx = [",I3," ",I3, "]; domy = [",I3," ",I3, "]; icell = [",F10.3," ",F10.3,"] jcell = [",F10.3," ",F10.3,"]")' & + , id, domx(1),domx(2), domy(1), domy(2), icell(domy(1)), icell(domy(2)-1), jcell(domx(1)), jcell(domx(2)-1) + print '("id = ", I2, " ids = [", I2, " ", I2, " ", I2, " ", I2,"]")', id, ids(1), ids(2), ids(3), ids(4) + ! libera a memória utilizada pelos processos não-raiz + ! pois eles não vão imprimir x e v pra csv + if (id /= 0) then + deallocate(v,x) + end if + ! aloca alguns vetores utilizados para o mpi_gatherv + allocate(rcounts(np), displs(np)) + + ! print '("Mesh divided in ", I2, "x", I2, " subregions." )', subx, suby + ! print*, 'ID ', id, 'iterv ', interv + call clean_mesh(malha, mesh, domx, domy,id,.false.) + ! print*, "bbb", id + if (id == 0) then + print*,'Press Return to continue' + ! read(*,*) + end if + call MPI_barrier(MPI_COMM_WORLD, ierr) + + ! id = 0 é o processo raiz + if (id == 0) call cpu_time(start) + + ! Agora dividimos para cada id uma região. + i = 0 + j = 1 + j2 = 1 + do while (t_fim > t) + ! print*, "L 1990" + call comp_F(GField, mesh,malha,propriedade,rcut,fric_term,domx,domy,ids,id,t) !altera Força + + ! IDS são os ids das regiões vizinhas, vetor(4) int posições, + ! DOMX e DOMY são vetores(2) int com o domínio (quat. de celulas) + ! que o processo vai cuidar (subdivisões) + ! print*, "L 1702", id + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, "L 1999", id + ! if (id == 0) print*, "TEMPO", t, '<<<<<<<<<<', id, i + call comp_x(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,domy,id, np) ! altera posição + ! print*, "L 2002", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + + call walls(icell,jcell,mesh,malha,domx,domy,wall(1:1),wall(2:2),wall(3:3),wall(4:4),id) ! altera posição e malha + + ! print*, "L 1714", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + call comp_F(GField, mesh,malha,propriedade,rcut,fric_term,domx,domy,ids,id,t) !altera força + ! print*, "L 1717" + call comp_v(malha,mesh,dt,t,propriedade,domx,domy) !altera velocidade + ! print*, "1721" + + t = t + dt + cont2 = cont2+1 + + if (i == interv_Td(j2)) then + call corr_K(malha,domx,domy,N,Td,propriedade, np,id,t) + j2 = j2 + 1 + end if + + if (i == interv(j)) then + deallocate(LT%lstrdb_N, LT%lstrdb_S, LT%lstrdb_E, LT%lstrdb_W, LT%lstrdb_D, & + LT%lstrint_N, LT%lstrint_S, LT%lstrint_E, LT%lstrint_W, LT%lstrint_D) + allocate(nxv(N*5),nxv_send(N*5)) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, 'MANDANDO PRINTAR', id, i + + call linked2vec(malha,domx,domy,nxv_send,aux1) + ! ! ! print*, "L nxv_send", nxv_send + !! print*, "LINKED 2 VEC, id", id, "aux1", aux1 + if (np > 1) then ! paralelo + ! call MPI_GATHER(sbuf, scount, MPI_integer, rbuf, rcount, MPI_integer, root, MPI_COMM_WORLD, ierr) + call MPI_GATHER(aux1, 1, MPI_integer, rcounts, 1, MPI_integer, 0, MPI_COMM_WORLD, ierr) + + ! print*, "Aux1 =", aux1, "id =", id + ! if (id == 0) print*, "rcounts =", rcounts + ! !print*, 'L IMPRE >', id + displs = 0 + + if (id == 0) then + do ii = 2,np + displs(ii) = displs(ii-1) + rcounts(ii-1) + end do + ! print*, "displs", displs + end if + ! print*,'E_tot0 = ',(comp_pot(mesh,malha,propriedade,rcut) +comp_K(malha,mesh,propriedade)) + ! MPI_GATHERV (sbuf, scount, stype, rbuf, rcounts, displs, rtype, root, comm, ierr) + call MPI_GATHERV(nxv_send, aux1, MPI_DOUBLE_PRECISION, nxv, rcounts, displs, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + ! o processo raiz imprime os resultados + ! print*, "id", id, "nxv_send:", nxv_send(1:aux1) + + ! if (id == 0) print*, "nxv:", nxv + + else + nxv = nxv_send + end if + + if (id == 0) then + do ii = 0, N-1 + + x(int(nxv(ii*5+1)),:) = [nxv(ii*5+2),nxv(ii*5+3)] + ! print*, 'x= ',[nxv(ii*5+2),nxv(ii*5+3)], 'i=',int(nxv(ii*5+1)) + v(int(nxv(ii*5+1)),:) = [nxv(ii*5+4),nxv(ii*5+5)] + end do + + call vec2csv(x,N,2,'position',j,t,nimpre,start) + call vec2csv(v,N,2,'velocity',j,t,nimpre,start) + + end if + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + j = j+1 + deallocate(nxv,nxv_send) + allocate(LT%lstrdb_N(NMPT*6),LT%lstrdb_S(NMPT*6),LT%lstrdb_E(NMPT*6),LT%lstrdb_W(NMPT*6),LT%lstrdb_D(NMPT*6), & + LT%lstrint_N(NMPT*4),LT%lstrint_S(NMPT*4),LT%lstrint_E(NMPT*4),LT%lstrint_W(NMPT*4), LT%lstrint_D(NMPT*4)) + end if + i = i+1 + ! print*, "L 1754 >", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end do + + ! print*, 'L 2112', id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + call clean_mesh(malha, mesh, domx, domy,id,.true.) + deallocate(LT%lstrdb_N, LT%lstrdb_S, LT%lstrdb_E, LT%lstrdb_W, LT%lstrdb_D, & + LT%lstrint_N, LT%lstrint_S, LT%lstrint_E, LT%lstrint_W, LT%lstrint_D) + if (id == 0) then + call cpu_time(finish) + 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." + close(22) + !call system('python csv2vtk_particles.py') + end if + + ! print*, 'L 1808' + call MPI_Finalize ( ierr ) +end program main \ No newline at end of file diff --git a/linkedlist.f90 b/linkedlist.f90 new file mode 100644 index 0000000..4c033af --- /dev/null +++ b/linkedlist.f90 @@ -0,0 +1,162 @@ +module linkedlist + implicit none + + private + + public :: list_t + public :: list_data + public :: list_init + public :: list_free + public :: list_remove + public :: list_insert + public :: list_put + public :: list_get + public :: list_next + public :: list_change + + ! A public variable to use as a MOLD for transfer() + integer, dimension(:), allocatable :: list_data + + ! Linked list node data type + type :: list_t + private + integer, dimension(:), pointer :: data => null() + type(list_t), pointer :: next => null() + end type list_t + +contains + + ! Initialize a head node SELF and optionally store the provided DATA. + subroutine list_init(self, data) + type(list_t), pointer :: self + integer, dimension(:), intent(in), optional :: data + + allocate(self) + nullify(self%next) + + if (present(data)) then + allocate(self%data(size(data))) + self%data = data + else + nullify(self%data) + end if + end subroutine list_init + + ! Free the entire list and all data, beginning at SELF + subroutine list_free(self) + type(list_t), pointer :: self + type(list_t), pointer :: current + type(list_t), pointer :: next + + current => self + do while (associated(current)) + next => current%next + if (associated(current%data)) then + deallocate(current%data) + nullify(current%data) + end if + deallocate(current) + nullify(current) + current => next + end do + + end subroutine list_free + + ! Remove element AFTER SELF + subroutine list_remove(self) + type(list_t), pointer :: self + ! type(list_t), pointer :: current + type(list_t), pointer :: next + + + + next => self%next + + self%next => next%next + ! print*, "A" + if (associated(next%data)) then + ! print*, "B" + deallocate(next%data) + ! print*, "C" + nullify(next%data) + ! print*, "D" + end if + + ! print*, "E" + deallocate(next) + ! print*, "F" + nullify(next) + end subroutine list_remove + + + ! Return the next node after SELF + function list_next(self) result(next) + type(list_t), pointer :: self + type(list_t), pointer :: next + next => self%next + end function list_next + + ! Insert a list node after SELF containing DATA (optional) + subroutine list_insert(self, data) + type(list_t), pointer :: self + integer, dimension(:), intent(in), optional :: data + type(list_t), pointer :: next + + allocate(next) + + if (present(data)) then + allocate(next%data(size(data))) + next%data = data + else + nullify(next%data) + end if + + next%next => self%next + self%next => next + + end subroutine list_insert + +! Change the NEXT node after a element CURRENT to after anOTHER (in other list for example). + subroutine list_change(current,other) + type(list_t), pointer :: current + type(list_t), pointer :: next + type(list_t), pointer :: other + + next => current%next + if (associated(next%next)) then + current%next => next%next + else + current%next => null() + end if + next%next => other%next + other%next => next + + + end subroutine list_change + + ! Store the encoded DATA in list node SELF + subroutine list_put(self, data) + type(list_t), pointer :: self + integer, dimension(:), intent(in) :: data + + if (associated(self%data)) then + deallocate(self%data) + nullify(self%data) + end if + self%data = data + end subroutine list_put + + ! Return the DATA stored in the node SELF + function list_get(self) result(data) + type(list_t), pointer :: self + integer, dimension(:), pointer :: data + data => self%data + end function list_get + +! function list_isnull(self) result(a) +! type(list_t), pointer :: self = .false. +! logical :: a +! if + + +end module linkedlist diff --git a/linkedtest.f90 b/linkedtest.f90 new file mode 100644 index 0000000..e3ec6fd --- /dev/null +++ b/linkedtest.f90 @@ -0,0 +1,134 @@ +! A simple generic linked list test program + +module data + implicit none + + private + public :: data_t + public :: data_ptr + + ! Data is stored in data_t + type :: data_t + real :: x + end type data_t + + ! A trick to allow us to store pointers in the list + type :: data_ptr + type(data_t), pointer :: p + end type data_ptr +end module data + +program list_test + use linkedlist + use data + implicit none + + integer :: cont + type(list_t), pointer :: ll => null(), node + type(data_t), target :: dat_a + type(data_t), target :: dat_b + type(data_ptr) :: ptr + + + print*, "enter" + read(*,*) + + do cont = 1, 100000 + ! Initialize two data objects + dat_a%x = 17.5 + dat_b%x = 3.0 + + ! Initialize the list with dat_a + allocate(ptr%p) + ptr%p%x = 17.5 ! dat_a%x + call list_init(ll, DATA=transfer(ptr, list_data)) + print *, 'Initializing list with data:', ptr%p + ! deallocate(ptr%p) + + allocate(ptr%p) + ! Insert dat_b into the list + ptr%p%x = 3.0 ! dat_b%x + call list_insert(ll, DATA=transfer(ptr, list_data)) + print *, 'Inserting node with data:', ptr%p + ! deallocate(ptr%p) + + allocate(ptr%p) + ! Insert dat_b into the list + ptr%p%x = 40.0 ! dat_b%x + call list_insert(ll, DATA=transfer(ptr, list_data)) + print *, 'Inserting node with data:', ptr%p + ! deallocate(ptr%p) + + dat_a%x = 0 + dat_b%x = 0 + + ! Get the head node + ptr = transfer(list_get(ll), ptr) + print *, 'Head node data:', ptr%p + deallocate(ptr%p) + ! Get the next node + node => list_next(ll) + ptr = transfer(list_get(list_next(ll)), ptr) + print *, 'Second node data:', ptr%p + deallocate(ptr%p) + ! Get the next node + node => list_next(node) + ptr = transfer(list_get(node), ptr) + print *, 'Third node data:', ptr%p + deallocate(ptr%p) + ! Free the list + + ! print*, "blablabla" + ! ! Get the head node + ! ptr = transfer(list_get(ll), ptr) + ! print *, 'Head node data:', ptr%p + ! deallocate(ptr%p) + ! ! Get the next node + ! node => list_next(ll) + ! ptr = transfer(list_get(list_next(ll)), ptr) + ! print *, 'Second node data:', ptr%p + ! deallocate(ptr%p) + ! ! Get the next node + ! node => list_next(node) + ! ptr = transfer(list_get(node), ptr) + ! print *, 'Third node data:', ptr%p + ! deallocate(ptr%p) + ! ! Free the list + + call list_free(ll) + end do + print*, "enter" + read(*,*) + +! print*, "enter" +! read(*,*) +! do cont = 1, 100000 +! ! Initialize two data objects +! dat_a%x = 17.5 +! dat_b%x = 3.0 + +! ! Initialize the list with dat_a +! ptr%p => dat_a +! call list_init(ll, DATA=transfer(ptr, list_data)) +! print *, 'Initializing list with data:', ptr%p + +! ! Insert dat_b into the list +! ptr%p => dat_b +! call list_insert(ll, DATA=transfer(ptr, list_data)) +! print *, 'Inserting node with data:', ptr%p + +! ! Get the head node +! ptr = transfer(list_get(ll), ptr) +! print *, 'Head node data:', ptr%p + +! ! Get the next node +! ptr = transfer(list_get(list_next(ll)), ptr) +! print *, 'Second node data:', ptr%p + +! ! Free the list +! call list_free(ll) +! end do + + + +end program list_test \ No newline at end of file diff --git a/m_config.f90 b/m_config.f90 new file mode 100644 index 0000000..05b0e7f --- /dev/null +++ b/m_config.f90 @@ -0,0 +1,1243 @@ +!> Module that allows working with a configuration file +module m_config + + implicit none + private + + !> The double precision kind-parameter + integer, parameter :: dp = kind(0.0d0) + + integer, parameter :: CFG_num_types = 4 !< Number of variable types + integer, parameter :: CFG_integer_type = 1 !< Integer type + integer, parameter :: CFG_real_type = 2 !< Real number type + integer, parameter :: CFG_string_type = 3 !< String type + integer, parameter :: CFG_logic_type = 4 !< Boolean/logical type + integer, parameter :: CFG_unknown_type = 0 !< Used before a variable is created + + !> Names of the types + character(len=10), parameter :: CFG_type_names(0:CFG_num_types) = & + [character(len=10) :: "storage", "integer", "real", "string ", "logical"] + + integer, parameter :: CFG_name_len = 80 !< Maximum length of variable names + integer, parameter :: CFG_string_len = 200 !< Fixed length of string type + + !> Maximum number of entries in a variable (if it's an array) + integer, parameter :: CFG_max_array_size = 1000 + + !> The separator(s) for array-like variables (space, comma, ', ", and tab) + character(len=*), parameter :: CFG_separators = " ,'"""//char(9) + + !> The separator for categories (stored in var_name) + character(len=*), parameter :: CFG_category_separator = "%" + + !> The default string for data that is not yet stored + character(len=*), parameter :: unstored_data_string = "__UNSTORED_DATA_STRING" + + !> The type of a configuration variable + type CFG_var_t + private + !> Name of the variable + character(len=CFG_name_len) :: var_name + !> Description of variable + character(len=CFG_string_len) :: description + !> Type of variable + integer :: var_type + !> Size of variable, 1 means scalar, > 1 means array + integer :: var_size + !> Whether the variable size is flexible + logical :: dynamic_size + !> Whether the variable's value has been requested + logical :: used + !> Data that has been read in for this variable + character(len=CFG_string_len) :: stored_data + + ! These are the arrays used for storage. In the future, a "pointer" based + ! approach could be used. + real(dp), allocatable :: real_data(:) + integer, allocatable :: int_data(:) + character(len=CFG_string_len), allocatable :: char_data(:) + logical, allocatable :: logic_data(:) + end type CFG_var_t + + !> The configuration that contains all the variables + type CFG_t + logical :: sorted = .false. + integer :: num_vars = 0 + type(CFG_var_t), allocatable :: vars(:) + end type CFG_t + + !> Interface to add variables to the configuration + interface CFG_add + module procedure add_real, add_real_array + module procedure add_int, add_int_array + module procedure add_string, add_string_array + module procedure add_logic, add_logic_array + end interface CFG_add + + !> Interface to get variables from the configuration + interface CFG_get + module procedure get_real, get_real_array + module procedure get_int, get_int_array + module procedure get_logic, get_logic_array + module procedure get_string, get_string_array + end interface CFG_get + + !> Interface to get variables from the configuration + interface CFG_add_get + module procedure add_get_real, add_get_real_array + module procedure add_get_int, add_get_int_array + module procedure add_get_logic, add_get_logic_array + module procedure add_get_string, add_get_string_array + end interface CFG_add_get + + ! Public types + public :: CFG_t + public :: CFG_integer_type + public :: CFG_real_type + public :: CFG_string_type + public :: CFG_logic_type + public :: CFG_type_names + + ! Constants + public :: CFG_name_len + public :: CFG_string_len + public :: CFG_max_array_size + + ! Public methods + public :: CFG_add + public :: CFG_get + public :: CFG_add_get + public :: CFG_get_size + public :: CFG_get_type + public :: CFG_check + public :: CFG_sort + public :: CFG_write + public :: CFG_write_markdown + public :: CFG_read_file + public :: CFG_update_from_arguments + public :: CFG_clear + +contains + + !> Read command line arguments. Both files and variables can be specified, for + !> example as: ./my_program config.cfg -n_runs=3 + subroutine CFG_update_from_arguments(cfg) + type(CFG_t),intent(inout) :: cfg + character(len=CFG_string_len) :: arg + integer :: ix + logical :: valid_syntax + + do ix = 1, command_argument_count() + call get_command_argument(ix, arg) + + if (arg(1:1) == '-') then + ! This sets a variable + call parse_line(cfg, arg(2:), valid_syntax) + if (.not. valid_syntax) then + call handle_error("Invalid variable specified on command line") + end if + else + ! This is a configuration files + call CFG_read_file(cfg, trim(arg)) + end if + end do + end subroutine CFG_update_from_arguments + + !> This routine will be called if an error occurs in one of the subroutines of + !> this module. + subroutine handle_error(err_string) + character(len=*), intent(in) :: err_string + + print *, "The following error occured in m_config:" + print *, trim(err_string) + + ! It is usually best to quit after an error, to make sure the error message + ! is not overlooked in the program's output + error stop + end subroutine handle_error + + !> Return the index of the variable with name 'var_name', or -1 if not found. + subroutine get_var_index(cfg, var_name, ix) + type(CFG_t), intent(in) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(out) :: ix + integer :: i + + if (cfg%sorted) then + call binary_search_variable(cfg, var_name, ix) + else + ! Linear search + do i = 1, cfg%num_vars + if (cfg%vars(i)%var_name == var_name) exit + end do + + ! If not found, set i to -1 + if (i == cfg%num_vars + 1) i = -1 + ix = i + end if + + end subroutine get_var_index + + !> Update the variables in the configartion with the values found in 'filename' + subroutine CFG_read_file(cfg, filename) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: filename + + integer, parameter :: my_unit = 123 + integer :: io_state + integer :: line_number + logical :: valid_syntax + character(len=CFG_name_len) :: line_fmt + character(len=CFG_string_len) :: err_string + character(len=CFG_string_len) :: line + character(len=CFG_name_len) :: category + + open(my_unit, file=trim(filename), status="old", action="read") + write(line_fmt, "(A,I0,A)") "(A", CFG_string_len, ")" + + category = "" ! Default category is empty + line_number = 0 + + do + read(my_unit, FMT=trim(line_fmt), ERR=998, end=999) line + line_number = line_number + 1 + + call parse_line(cfg, line, valid_syntax, category) + + if (.not. valid_syntax) then + write(err_string, *) "Cannot read line ", line_number, & + " from ", trim(filename) + call handle_error(err_string) + end if + end do + + ! Error handling +998 write(err_string, "(A,I0,A,I0)") " IOSTAT = ", io_state, & + " while reading from " // trim(filename) // " at line ", & + line_number + call handle_error("CFG_read_file:" // err_string) + + ! Routine ends here if the end of "filename" is reached +999 close(my_unit, iostat=io_state) + + end subroutine CFG_read_file + + !> Update the cfg by parsing one line + subroutine parse_line(cfg, line_arg, valid_syntax, category_arg) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: line_arg + logical, intent(out) :: valid_syntax + character(len=CFG_name_len), intent(inout), optional :: category_arg + character(len=CFG_name_len) :: var_name, category + integer :: ix, equal_sign_ix + logical :: append + character(len=CFG_string_len) :: line + + valid_syntax = .true. + + ! Work on a copy + line = line_arg + category = "" + if (present(category_arg)) category = category_arg + + call trim_comment(line, '#;') + + ! Skip empty lines + if (line == "") return + + ! Locate the '=' sign + equal_sign_ix = scan(line, '=') + + ! if there is no '='-sign then a category is indicated + if (equal_sign_ix == 0) then + line = adjustl(line) + + ! The category name should appear like this: [category_name] + ix = scan(line, ']') + if (line(1:1) /= '[' .or. ix == 0) then + valid_syntax = .false. + return + else + if (present(category_arg)) category_arg = line(2:ix-1) + return + end if + end if + + if (line(equal_sign_ix-1:equal_sign_ix) == '+=') then + append = .true. + var_name = line(1 : equal_sign_ix - 2) ! Set variable name + else + append = .false. + var_name = line(1 : equal_sign_ix - 1) ! Set variable name + end if + + ! If there are less than two spaces or a tab, reset to no category + if (var_name(1:2) /= " " .and. var_name(1:1) /= char(9)) then + category = "" + end if + + ! Remove leading blanks + var_name = adjustl(var_name) + + ! Add category if it is defined + if (category /= "") then + var_name = trim(category) // CFG_category_separator // var_name + end if + + line = line(equal_sign_ix + 1:) ! Set line to the values behind the '=' sign + + ! Find variable corresponding to name in file + call get_var_index(cfg, var_name, ix) + + if (ix <= 0) then + ! Variable still needs to be created, for now store data as a string + call prepare_store_var(cfg, trim(var_name), CFG_unknown_type, 1, & + "Not yet created", ix, .false.) + cfg%vars(ix)%stored_data = line + else + if (append) then + cfg%vars(ix)%stored_data = & + trim(cfg%vars(ix)%stored_data) // trim(line) + else + cfg%vars(ix)%stored_data = line + end if + + ! If type is known, read in values + if (cfg%vars(ix)%var_type /= CFG_unknown_type) then + call read_variable(cfg%vars(ix)) + end if + end if + + end subroutine parse_line + + subroutine read_variable(var) + type(CFG_var_t), intent(inout) :: var + integer :: n, n_entries + integer :: ix_start(CFG_max_array_size) + integer :: ix_end(CFG_max_array_size), stat + + ! Get the start and end positions of the line content, and the number of entries + call get_fields_string(var%stored_data, CFG_separators, & + CFG_max_array_size, n_entries, ix_start, ix_end) + + if (var%var_size /= n_entries) then + + if (.not. var%dynamic_size) then + ! Allow strings of length 1 to be automatically concatenated + if (var%var_type == CFG_string_type .and. var%var_size == 1) then + var%char_data(1) = trim(var%stored_data(ix_start(1):ix_end(n_entries))) + return ! Leave routine + else + call handle_error("read_variable: variable [" // & + & trim(var%var_name) // "] has the wrong size") + end if + else + var%var_size = n_entries + call resize_storage(var) + end if + end if + + do n = 1, n_entries + stat = 0 + select case (var%var_type) + case (CFG_integer_type) + read(var%stored_data(ix_start(n):ix_end(n)), *, iostat=stat) var%int_data(n) + case (CFG_real_type) + read(var%stored_data(ix_start(n):ix_end(n)), *, iostat=stat) var%real_data(n) + case (CFG_string_type) + var%char_data(n) = trim(var%stored_data(ix_start(n):ix_end(n))) + case (CFG_logic_type) + read(var%stored_data(ix_start(n):ix_end(n)), *, iostat=stat) var%logic_data(n) + end select + + if(stat /= 0) then + write (*, *) "** m_config error **" + write (*, *) "reading variable: ", trim(var%var_name) + write (*, *) "variable type: ", trim(CFG_type_names(var%var_type)) + write (*, *) "parsing value: ", var%stored_data(ix_start(n):ix_end(n)) + write (*, "(A,I0)") " iostat value: ", stat + stop + endif + end do + end subroutine read_variable + + subroutine trim_comment(line, comment_chars) + character(len=*), intent(inout) :: line + character(len=*), intent(in) :: comment_chars + character :: current_char, need_char + integer :: n + + ! Strip comments, but only outside quoted strings (so that var = '#yolo' is + ! valid when # is a comment char) + need_char = "" + + do n = 1, len(line) + current_char = line(n:n) + + if (need_char == "") then + if (current_char == "'") then + need_char = "'" ! Open string + else if (current_char == '"') then + need_char = '"' ! Open string + else if (index(comment_chars, current_char) /= 0) then + line = line(1:n-1) ! Trim line up to comment character + exit + end if + else if (current_char == need_char) then + need_char = "" ! Close string + end if + + end do + + end subroutine trim_comment + + subroutine CFG_check(cfg) + type(CFG_t), intent(in) :: cfg + integer :: n + character(len=CFG_string_len) :: err_string + + do n = 1, cfg%num_vars + if (cfg%vars(n)%var_type == CFG_unknown_type) then + write(err_string, *) "CFG_check: unknown variable ", & + trim(cfg%vars(n)%var_name), " specified" + call handle_error(err_string) + end if + end do + end subroutine CFG_check + + !> This routine writes the current configuration to a file with descriptions + subroutine CFG_write(cfg_in, filename, hide_unused) + use iso_fortran_env + type(CFG_t), intent(in) :: cfg_in + character(len=*), intent(in) :: filename + logical, intent(in), optional :: hide_unused + logical :: hide_not_used + type(CFG_t) :: cfg + integer :: i, j, io_state, myUnit + character(len=CFG_name_len) :: name_format, var_name + character(len=CFG_name_len) :: category, prev_category + character(len=CFG_string_len) :: err_string + + hide_not_used = .false. + if (present(hide_unused)) hide_not_used = hide_unused + + ! Always print a sorted configuration + cfg = cfg_in + if (.not. cfg%sorted) call CFG_sort(cfg) + + write(name_format, FMT="(A,I0,A)") "(A,A", CFG_name_len, ",A)" + + if (filename == "stdout") then + myUnit = output_unit + else + myUnit = 333 + open(myUnit, FILE=filename, ACTION="WRITE") + end if + + category = "" + prev_category = "" + + do i = 1, cfg%num_vars + if (.not. cfg%vars(i)%used .and. hide_not_used) cycle + if (cfg%vars(i)%var_type == CFG_unknown_type) cycle + + ! Write category when it changes + call split_category(cfg%vars(i), category, var_name) + + if (category /= prev_category .and. category /= '') then + write(myUnit, ERR=998, FMT="(A)") '[' // trim(category) // ']' + prev_category = category + end if + + ! Indent if inside category + if (category /= "") then + write(myUnit, ERR=998, FMT="(A,A,A)") " # ", & + trim(cfg%vars(i)%description), ":" + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A)") & + " " // trim(var_name) // " =" + else + write(myUnit, ERR=998, FMT="(A,A,A)") "# ", & + trim(cfg%vars(i)%description), ":" + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A)") & + trim(var_name) // " =" + end if + + select case(cfg%vars(i)%var_type) + case (CFG_integer_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A,I0)") & + " ", cfg%vars(i)%int_data(j) + end do + case (CFG_real_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A,E11.4)") & + " ", cfg%vars(i)%real_data(j) + end do + case (CFG_string_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A)") & + " '" // trim(cfg%vars(i)%char_data(j)) // "'" + end do + case (CFG_logic_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A,L1)") & + " ", cfg%vars(i)%logic_data(j) + end do + end select + write(myUnit, ERR=998, FMT="(A)") "" + write(myUnit, ERR=998, FMT="(A)") "" + end do + + if (myUnit /= output_unit) close(myUnit, ERR=999, IOSTAT=io_state) + call CFG_check(cfg_in) + return + +998 continue + write(err_string, *) "CFG_write error: io_state = ", io_state, & + " while writing ", trim(var_name), " to ", filename + call handle_error(err_string) + +999 continue ! If there was an error, the routine will end here + write(err_string, *) "CFG_write error: io_state = ", io_state, & + " while writing to ", filename + call handle_error(err_string) + + end subroutine CFG_write + + !> This routine writes the current configuration to a markdown file + subroutine CFG_write_markdown(cfg_in, filename, hide_unused) + use iso_fortran_env + type(CFG_t), intent(in) :: cfg_in + character(len=*), intent(in) :: filename + logical, intent(in), optional :: hide_unused + logical :: hide_not_used + integer :: i, j, io_state, myUnit + type(CFG_t) :: cfg + character(len=CFG_name_len) :: name_format, var_name + character(len=CFG_name_len) :: category, prev_category + character(len=CFG_string_len) :: err_string + + hide_not_used = .false. + if (present(hide_unused)) hide_not_used = hide_unused + + ! Always print a sorted configuration + cfg = cfg_in + if (.not. cfg%sorted) call CFG_sort(cfg) + + write(name_format, FMT="(A,I0,A)") "(A,A", CFG_name_len, ",A)" + + if (filename == "stdout") then + myUnit = output_unit + else + myUnit = 333 + open(myUnit, FILE=filename, ACTION="WRITE") + end if + + category = "" + prev_category = "X" + write(myUnit, ERR=998, FMT="(A)") "# Configuration file (markdown format)" + write(myUnit, ERR=998, FMT="(A)") "" + + do i = 1, cfg%num_vars + + if (.not. cfg%vars(i)%used .and. hide_not_used) cycle + if (cfg%vars(i)%var_type == CFG_unknown_type) cycle + + ! Write category when it changes + call split_category(cfg%vars(i), category, var_name) + + if (category /= prev_category) then + if (category == "") category = "No category" + write(myUnit, ERR=998, FMT="(A)") '## ' // trim(category) + write(myUnit, ERR=998, FMT="(A)") "" + prev_category = category + end if + + write(myUnit, ERR=998, FMT="(A)") "* " // trim(cfg%vars(i)%description) + write(myUnit, ERR=998, FMT="(A)") "" + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A)") & + ' ' // trim(var_name) // " =" + + select case(cfg%vars(i)%var_type) + case (CFG_integer_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A,I0)") & + " ", cfg%vars(i)%int_data(j) + end do + case (CFG_real_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A,E11.4)") & + " ", cfg%vars(i)%real_data(j) + end do + case (CFG_string_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A)") & + " '" // trim(cfg%vars(i)%char_data(j)) // "'" + end do + case (CFG_logic_type) + do j = 1, cfg%vars(i)%var_size + write(myUnit, ADVANCE="NO", ERR=998, FMT="(A,L1)") & + " ", cfg%vars(i)%logic_data(j) + end do + end select + write(myUnit, ERR=998, FMT="(A)") "" + write(myUnit, ERR=998, FMT="(A)") "" + end do + + if (myUnit /= output_unit) close(myUnit, ERR=999, IOSTAT=io_state) + call CFG_check(cfg_in) + return + +998 continue + write(err_string, *) "CFG_write_markdown error: io_state = ", io_state, & + " while writing ", trim(var_name), " to ", filename + call handle_error(err_string) + +999 continue ! If there was an error, the routine will end here + write(err_string, *) "CFG_write_markdown error: io_state = ", io_state, & + " while writing to ", filename + call handle_error(err_string) + + end subroutine CFG_write_markdown + + subroutine split_category(variable, category, var_name) + type(CFG_var_t), intent(in) :: variable + character(CFG_name_len), intent(out) :: category + character(CFG_name_len), intent(out) :: var_name + integer :: ix + + ix = index(variable%var_name, CFG_category_separator) + + if (ix == 0) then + category = "" + var_name = variable%var_name + else + category = variable%var_name(1:ix-1) + var_name = variable%var_name(ix+1:) + end if + + end subroutine split_category + + !> Resize the storage size of variable, which can be of type integer, logical, + !> real or character + subroutine resize_storage(variable) + type(CFG_var_t), intent(inout) :: variable + + select case (variable%var_type) + case (CFG_integer_type) + deallocate( variable%int_data ) + allocate( variable%int_data(variable%var_size) ) + case (CFG_logic_type) + deallocate( variable%logic_data ) + allocate( variable%logic_data(variable%var_size) ) + case (CFG_real_type) + deallocate( variable%real_data ) + allocate( variable%real_data(variable%var_size) ) + case (CFG_string_type) + deallocate( variable%char_data ) + allocate( variable%char_data(variable%var_size) ) + end select + end subroutine resize_storage + + !> Helper routine to store variables. This is useful because a lot of the same + !> code is executed for the different types of variables. + subroutine prepare_store_var(cfg, var_name, var_type, var_size, & + description, ix, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, description + integer, intent(in) :: var_type, var_size + integer, intent(out) :: ix !< Index of variable + logical, intent(in), optional :: dynamic_size + + ! Check if variable already exists + call get_var_index(cfg, var_name, ix) + + if (ix == -1) then ! Create a new variable + call ensure_free_storage(cfg) + cfg%sorted = .false. + ix = cfg%num_vars + 1 + cfg%num_vars = cfg%num_vars + 1 + cfg%vars(ix)%used = .false. + cfg%vars(ix)%stored_data = unstored_data_string + else + ! Only allowed when the variable is not yet created + if (cfg%vars(ix)%var_type /= CFG_unknown_type) then + call handle_error("prepare_store_var: variable [" // & + & trim(var_name) // "] already exists") + end if + end if + + cfg%vars(ix)%var_name = var_name + cfg%vars(ix)%description = description + cfg%vars(ix)%var_type = var_type + cfg%vars(ix)%var_size = var_size + + if (present(dynamic_size)) then + cfg%vars(ix)%dynamic_size = dynamic_size + else + cfg%vars(ix)%dynamic_size = .false. + end if + + select case (var_type) + case (CFG_integer_type) + allocate( cfg%vars(ix)%int_data(var_size) ) + case (CFG_real_type) + allocate( cfg%vars(ix)%real_data(var_size) ) + case (CFG_string_type) + allocate( cfg%vars(ix)%char_data(var_size) ) + case (CFG_logic_type) + allocate( cfg%vars(ix)%logic_data(var_size) ) + end select + + end subroutine prepare_store_var + + !> Helper routine to get variables. This is useful because a lot of the same + !> code is executed for the different types of variables. + subroutine prepare_get_var(cfg, var_name, var_type, var_size, ix) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(in) :: var_type, var_size + integer, intent(out) :: ix + character(len=CFG_string_len) :: err_string + + call get_var_index(cfg, var_name, ix) + + if (ix == -1) then + call handle_error("CFG_get: variable ["//var_name//"] not found") + else if (cfg%vars(ix)%var_type /= var_type) then + write(err_string, fmt="(A)") "CFG_get: variable [" & + // var_name // "] has different type (" // & + trim(CFG_type_names(cfg%vars(ix)%var_type)) // & + ") than requested (" // trim(CFG_type_names(var_type)) // ")" + call handle_error(err_string) + else if (cfg%vars(ix)%var_size /= var_size) then + write(err_string, fmt="(A,I0,A,I0,A)") "CFG_get: variable [" & + // var_name // "] has different size (", cfg%vars(ix)%var_size, & + ") than requested (", var_size, ")" + call handle_error(err_string) + else ! All good, variable will be used + cfg%vars(ix)%used = .true. + end if + end subroutine prepare_get_var + + !> Add a configuration variable with a real value + subroutine add_real(cfg, var_name, real_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + real(dp), intent(in) :: real_data + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_real_type, 1, comment, ix) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%real_data(1) = real_data + end if + end subroutine add_real + + !> Add a configuration variable with an array of type + ! real + subroutine add_real_array(cfg, var_name, real_data, comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + real(dp), intent(in) :: real_data(:) + logical, intent(in), optional :: dynamic_size + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_real_type, & + size(real_data), comment, ix, dynamic_size) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%real_data = real_data + end if + end subroutine add_real_array + + !> Add a configuration variable with an integer value + subroutine add_int(cfg, var_name, int_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + integer, intent(in) :: int_data + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_integer_type, 1, comment, ix) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%int_data(1) = int_data + end if + end subroutine add_int + + !> Add a configuration variable with an array of type integer + subroutine add_int_array(cfg, var_name, int_data, comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + integer, intent(in) :: int_data(:) + logical, intent(in), optional :: dynamic_size + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_integer_type, & + size(int_data), comment, ix, dynamic_size) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%int_data = int_data + end if + end subroutine add_int_array + + !> Add a configuration variable with an character value + subroutine add_string(cfg, var_name, char_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment, char_data + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_string_type, 1, comment, ix) + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%char_data(1) = char_data + end if + end subroutine add_string + + !> Add a configuration variable with an array of type character + subroutine add_string_array(cfg, var_name, char_data, & + comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment, char_data(:) + logical, intent(in), optional :: dynamic_size + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_string_type, & + size(char_data), comment, ix, dynamic_size) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%char_data = char_data + end if + end subroutine add_string_array + + !> Add a configuration variable with an logical value + subroutine add_logic(cfg, var_name, logic_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + logical, intent(in) :: logic_data + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_logic_type, 1, comment, ix) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%logic_data(1) = logic_data + end if + end subroutine add_logic + + !> Add a configuration variable with an array of type logical + subroutine add_logic_array(cfg, var_name, logic_data, & + comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + logical, intent(in) :: logic_data(:) + logical, intent(in), optional :: dynamic_size + integer :: ix + + call prepare_store_var(cfg, var_name, CFG_logic_type, & + size(logic_data), comment, ix, dynamic_size) + + if (cfg%vars(ix)%stored_data /= unstored_data_string) then + call read_variable(cfg%vars(ix)) + else + cfg%vars(ix)%logic_data = logic_data + end if + end subroutine add_logic_array + + !> Get a real array of a given name + subroutine get_real_array(cfg, var_name, real_data) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + real(dp), intent(inout) :: real_data(:) + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_real_type, & + size(real_data), ix) + real_data = cfg%vars(ix)%real_data + end subroutine get_real_array + + !> Get a integer array of a given name + subroutine get_int_array(cfg, var_name, int_data) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(inout) :: int_data(:) + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_integer_type, & + size(int_data), ix) + int_data = cfg%vars(ix)%int_data + end subroutine get_int_array + + !> Get a character array of a given name + subroutine get_string_array(cfg, var_name, char_data) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + character(len=*), intent(inout) :: char_data(:) + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_string_type, & + size(char_data), ix) + char_data = cfg%vars(ix)%char_data + end subroutine get_string_array + + !> Get a logical array of a given name + subroutine get_logic_array(cfg, var_name, logic_data) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + logical, intent(inout) :: logic_data(:) + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_logic_type, & + size(logic_data), ix) + logic_data = cfg%vars(ix)%logic_data + end subroutine get_logic_array + + !> Get a real value of a given name + subroutine get_real(cfg, var_name, res) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + real(dp), intent(out) :: res + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_real_type, 1, ix) + res = cfg%vars(ix)%real_data(1) + end subroutine get_real + + !> Get a integer value of a given name + subroutine get_int(cfg, var_name, res) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(inout) :: res + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_integer_type, 1, ix) + res = cfg%vars(ix)%int_data(1) + end subroutine get_int + + !> Get a logical value of a given name + subroutine get_logic(cfg, var_name, res) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + logical, intent(out) :: res + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_logic_type, 1, ix) + res = cfg%vars(ix)%logic_data(1) + end subroutine get_logic + + !> Get a character value of a given name + subroutine get_string(cfg, var_name, res) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name + character(len=*), intent(out) :: res + integer :: ix + + call prepare_get_var(cfg, var_name, CFG_string_type, 1, ix) + res = cfg%vars(ix)%char_data(1) + end subroutine get_string + + !> Get or add a real array of a given name + subroutine add_get_real_array(cfg, var_name, real_data, & + comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + real(dp), intent(inout) :: real_data(:) + logical, intent(in), optional :: dynamic_size + + call add_real_array(cfg, var_name, real_data, comment, dynamic_size) + call get_real_array(cfg, var_name, real_data) + end subroutine add_get_real_array + + !> Get or add a integer array of a given name + subroutine add_get_int_array(cfg, var_name, int_data, & + comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + integer, intent(inout) :: int_data(:) + logical, intent(in), optional :: dynamic_size + + call add_int_array(cfg, var_name, int_data, comment, dynamic_size) + call get_int_array(cfg, var_name, int_data) + end subroutine add_get_int_array + + !> Get or add a character array of a given name + subroutine add_get_string_array(cfg, var_name, char_data, & + comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + character(len=*), intent(inout) :: char_data(:) + logical, intent(in), optional :: dynamic_size + + call add_string_array(cfg, var_name, char_data, comment, dynamic_size) + call get_string_array(cfg, var_name, char_data) + end subroutine add_get_string_array + + !> Get or add a logical array of a given name + subroutine add_get_logic_array(cfg, var_name, logic_data, & + comment, dynamic_size) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + logical, intent(inout) :: logic_data(:) + logical, intent(in), optional :: dynamic_size + + call add_logic_array(cfg, var_name, logic_data, comment, dynamic_size) + call get_logic_array(cfg, var_name, logic_data) + end subroutine add_get_logic_array + + !> Get or add a real value of a given name + subroutine add_get_real(cfg, var_name, real_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + real(dp), intent(inout) :: real_data + + call add_real(cfg, var_name, real_data, comment) + call get_real(cfg, var_name, real_data) + end subroutine add_get_real + + !> Get or add a integer value of a given name + subroutine add_get_int(cfg, var_name, int_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + integer, intent(inout) :: int_data + + call add_int(cfg, var_name, int_data, comment) + call get_int(cfg, var_name, int_data) + end subroutine add_get_int + + !> Get or add a logical value of a given name + subroutine add_get_logic(cfg, var_name, logical_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + logical, intent(inout) :: logical_data + + call add_logic(cfg, var_name, logical_data, comment) + call get_logic(cfg, var_name, logical_data) + end subroutine add_get_logic + + !> Get a character value of a given name + subroutine add_get_string(cfg, var_name, string_data, comment) + type(CFG_t), intent(inout) :: cfg + character(len=*), intent(in) :: var_name, comment + character(len=*), intent(inout) :: string_data + + call add_string(cfg, var_name, string_data, comment) + call get_string(cfg, var_name, string_data) + end subroutine add_get_string + + !> Get the size of a variable + subroutine CFG_get_size(cfg, var_name, res) + type(CFG_t), intent(in) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(out) :: res + integer :: ix + + call get_var_index(cfg, var_name, ix) + if (ix /= -1) then + res = cfg%vars(ix)%var_size + else + res = -1 + call handle_error("CFG_get_size: variable ["//var_name//"] not found") + end if + end subroutine CFG_get_size + + !> Get the type of a given variable of a configuration type + subroutine CFG_get_type(cfg, var_name, res) + type(CFG_t), intent(in) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(out) :: res + integer :: ix + + call get_var_index(cfg, var_name, ix) + + if (ix /= -1) then + res = cfg%vars(ix)%var_type + else + res = -1 + call handle_error("CFG_get_type: variable ["//var_name//"] not found") + end if + end subroutine CFG_get_type + + !> Routine to ensure that enough storage is allocated for the configuration + !> type. If not the new size will be twice as much as the current size. If no + !> storage is allocated yet a minumum amount of starage is allocated. + subroutine ensure_free_storage(cfg) + type(CFG_t), intent(inout) :: cfg + type(CFG_var_t), allocatable :: cfg_copy(:) + integer, parameter :: min_dyn_size = 100 + integer :: cur_size, new_size + + if (allocated(cfg%vars)) then + cur_size = size(cfg%vars) + + if (cur_size < cfg%num_vars + 1) then + new_size = 2 * cur_size + allocate(cfg_copy(cur_size)) + cfg_copy = cfg%vars + deallocate(cfg%vars) + allocate(cfg%vars(new_size)) + cfg%vars(1:cur_size) = cfg_copy + end if + else + allocate(cfg%vars(min_dyn_size)) + end if + + end subroutine ensure_free_storage + + !> Routine to find the indices of entries in a string + subroutine get_fields_string(line, delims, n_max, n_found, ixs_start, ixs_end) + !> The line from which we want to read + character(len=*), intent(in) :: line + !> A string with delimiters. For example delims = " ,'"""//char(9) + character(len=*), intent(in) :: delims + !> Maximum number of entries to read in + integer, intent(in) :: n_max + !> Number of entries found + integer, intent(inout) :: n_found + !> On return, ix_start(i) holds the starting point of entry i + integer, intent(inout) :: ixs_start(n_max) + !> On return, ix_end(i) holds the end point of entry i + integer, intent(inout) :: ixs_end(n_max) + + integer :: ix, ix_prev + + ix_prev = 0 + n_found = 0 + + do while (n_found < n_max) + + ! Find the starting point of the next entry (a non-delimiter value) + ix = verify(line(ix_prev+1:), delims) + if (ix == 0) exit + + n_found = n_found + 1 + ixs_start(n_found) = ix_prev + ix ! This is the absolute position in 'line' + + ! Get the end point of the current entry (next delimiter index minus one) + ix = scan(line(ixs_start(n_found)+1:), delims) - 1 + + if (ix == -1) then ! If there is no last delimiter, + ixs_end(n_found) = len(line) ! the end of the line is the endpoint + else + ixs_end(n_found) = ixs_start(n_found) + ix + end if + + ix_prev = ixs_end(n_found) ! We continue to search from here + end do + + end subroutine get_fields_string + + !> Performa a binary search for the variable 'var_name' + subroutine binary_search_variable(cfg, var_name, ix) + type(CFG_t), intent(in) :: cfg + character(len=*), intent(in) :: var_name + integer, intent(out) :: ix + integer :: i_min, i_max, i_mid + + i_min = 1 + i_max = cfg%num_vars + ix = - 1 + + do while (i_min < i_max) + i_mid = i_min + (i_max - i_min) / 2 + if ( llt(cfg%vars(i_mid)%var_name, var_name) ) then + i_min = i_mid + 1 + else + i_max = i_mid + end if + end do + + ! If not found, binary_search_variable is not set here, and stays -1 + if (i_max == i_min .and. cfg%vars(i_min)%var_name == var_name) then + ix = i_min + else + ix = -1 + end if + end subroutine binary_search_variable + + !> Sort the variables for faster lookup + subroutine CFG_sort(cfg) + type(CFG_t), intent(inout) :: cfg + + call qsort_config(cfg%vars(1:cfg%num_vars)) + cfg%sorted = .true. + end subroutine CFG_sort + + !> Simple implementation of quicksort algorithm to sort the variable list alphabetically. + recursive subroutine qsort_config(list) + type(CFG_var_t), intent(inout) :: list(:) + integer :: split_pos + + if (size(list) > 1) then + call parition_var_list(list, split_pos) + call qsort_config( list(:split_pos-1) ) + call qsort_config( list(split_pos:) ) + end if + end subroutine qsort_config + + !> Helper routine for quicksort, to perform partitioning + subroutine parition_var_list(list, marker) + type(CFG_var_t), intent(inout) :: list(:) + integer, intent(out) :: marker + integer :: left, right, pivot_ix + type(CFG_var_t) :: temp + character(len=CFG_name_len) :: pivot_value + + left = 0 + right = size(list) + 1 + + ! Take the middle element as pivot + pivot_ix = size(list) / 2 + pivot_value = list(pivot_ix)%var_name + + do while (left < right) + + right = right - 1 + do while (lgt(list(right)%var_name, pivot_value)) + right = right - 1 + end do + + left = left + 1 + do while (lgt(pivot_value, list(left)%var_name)) + left = left + 1 + end do + + if (left < right) then + temp = list(left) + list(left) = list(right) + list(right) = temp + end if + end do + + if (left == right) then + marker = left + 1 + else + marker = left + end if + end subroutine parition_var_list + + !> Clear all data from a CFG_t object, so that it can be reused. Note that + !> this also happens automatically when such an object goes out of scope. + subroutine CFG_clear(cfg) + implicit none + type(CFG_t) :: cfg + + cfg%sorted = .false. + cfg%num_vars = 0 + if(allocated(cfg%vars)) then + deallocate(cfg%vars) + endif + end subroutine CFG_clear + +end module m_config diff --git a/makefile b/makefile new file mode 100644 index 0000000..ea877c3 --- /dev/null +++ b/makefile @@ -0,0 +1,39 @@ +# Start of the makefile +# Defining variables +objects = m_config.o linkedlist.o lennard.o data.o matprint.o randnormal.o mod1.o +f90comp = mpiifort +switch = -O3 +# Makefile +execname: $(objects) + $(f90comp) -o lennard.out $(switch) $(objects) +randnormal.mod: randnormal.o randnormal.f90 + $(f90comp) -c $(switch) randnormal.f90 +mod1.mod: mod1.o mod1.f90 + $(f90comp) -c $(switch) mod1.f90 +matprint.mod: matprint.o matprint.f90 + $(f90comp) -c $(switch) matprint.f90 +data.mod: data.o data.f90 + $(f90comp) -c $(switch) data.f90 +linkedlist.mod: linkedlist.o linkedlist.f90 + $(f90comp) -c $(switch) linkedlist.f90 +m_config.mod: m_config.o m_config.f90 + $(f90comp) -c $(switch) m_config.f90 +randnormal.o: randnormal.f90 + $(f90comp) -c $(switch) randnormal.f90 +mod1.o: mod1.f90 + $(f90comp) -c $(switch) mod1.f90 +matprint.o: matprint.f90 + $(f90comp) -c $(switch) matprint.f90 +data.o: data.f90 + $(f90comp) -c $(switch) data.f90 +linkedlist.o: linkedlist.f90 + $(f90comp) -c $(switch) linkedlist.f90 +m_config.o: m_config.f90 + $(f90comp) -c $(switch) m_config.f90 +lennard.o: mod1.mod linkedlist.mod data.mod matprint.mod m_config.mod randnormal.mod lennard.f90 + $(f90comp) -c $(switch) lennard.f90 +# Cleaning everything +clean: + rm $(objects) + rm -f *.o *.mod *.MOD +# End of the makefile diff --git a/matprint.f90 b/matprint.f90 new file mode 100644 index 0000000..c1bb9bc --- /dev/null +++ b/matprint.f90 @@ -0,0 +1,26 @@ +! Este módulo imprime a matriz +module matprint + contains + subroutine printm(m,r) + use mod1 + real(dp), dimension(:,:), intent(in) :: m + integer,dimension(2) :: s + integer :: i,j + character(*), optional, intent(in) :: r + + s = shape(m) + if(present(r)) then + print*, r,' =' + else + print*, ' =' + end if + + do i = 1,s(1) + do j = 1,s(2) + write(*,'(EN14.3)',advance='no') m(i,j) + end do + write(*,*) + end do + write(*,*) ' ' + end subroutine printm +end module matprint diff --git a/mod1.f90 b/mod1.f90 new file mode 100644 index 0000000..a195a54 --- /dev/null +++ b/mod1.f90 @@ -0,0 +1,9 @@ +module mod1 +! use linkedlist + integer, parameter:: dp=kind(0.d0) ! double precision + ! vetor de strings + type string + character(len=:), allocatable :: str + integer :: quanty + end type string +end module mod1 diff --git a/normal.f90 b/normal.f90 new file mode 100644 index 0000000..012dd6d --- /dev/null +++ b/normal.f90 @@ -0,0 +1,1583 @@ +function c4_normal_01 ( seed ) + +!*****************************************************************************80 +! +!! C4_NORMAL_01 returns a unit pseudonormal C4. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, complex ( kind = 4 ) C4_NORMAL_01, a unit pseudonormal value. +! + implicit none + + complex ( kind = 4 ) c4_normal_01 + real ( kind = 4 ), parameter :: r4_pi = 3.141592653589793E+00 + real ( kind = 4 ) r4_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 4 ) v1 + real ( kind = 4 ) v2 + real ( kind = 4 ) x_c + real ( kind = 4 ) x_r + + v1 = r4_uniform_01 ( seed ) + v2 = r4_uniform_01 ( seed ) + + x_r = sqrt ( - 2.0E+00 * log ( v1 ) ) * cos ( 2.0E+00 * r4_pi * v2 ) + x_c = sqrt ( - 2.0E+00 * log ( v1 ) ) * sin ( 2.0E+00 * r4_pi * v2 ) + + c4_normal_01 = cmplx ( x_r, x_c ) + + return +end +function c8_normal_01 ( seed ) + +!*****************************************************************************80 +! +!! C8_NORMAL_01 returns a unit pseudonormal C8. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) SEED, a seed for the random number +! generator. +! +! Output, complex ( kind = 8 ) C8_NORMAL_01, a sample of the PDF. +! +! Output, integer ( kind = 4 ) SEED, a seed for the random number +! generator. +! + implicit none + + complex ( kind = 8 ) c8_normal_01 + real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00 + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 8 ) v1 + real ( kind = 8 ) v2 + real ( kind = 8 ) x_c + real ( kind = 8 ) x_r + + v1 = r8_uniform_01 ( seed ) + v2 = r8_uniform_01 ( seed ) + + x_r = sqrt ( - 2.0D+00 * log ( v1 ) ) * cos ( 2.0D+00 * r8_pi * v2 ) + x_c = sqrt ( - 2.0D+00 * log ( v1 ) ) * sin ( 2.0D+00 * r8_pi * v2 ) + + c8_normal_01 = cmplx ( x_r, x_c, kind = 8 ) + + return +end +function i4_normal_ab ( a, b, seed ) + +!*****************************************************************************80 +! +!! I4_NORMAL_AB returns a scaled pseudonormal I4. +! +! Discussion: +! +! The normal probability distribution function (PDF) is sampled, +! with mean A and standard deviation B. +! +! The result is then rounded to the nearest integer. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, real ( kind = 4 ) A, the mean of the PDF. +! +! Input, real ( kind = 4 ) B, the standard deviation of the PDF. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the +! random number generator. +! +! Output, integer ( kind = 4 ) I4_NORMAL_AB, a sample of the normal PDF. +! + implicit none + + real ( kind = 4 ) a + real ( kind = 4 ) b + integer ( kind = 4 ) i4_normal_ab + real ( kind = 4 ) r1 + real ( kind = 4 ) r2 + real ( kind = 4 ), parameter :: r4_pi = 3.141592653589793E+00 + real ( kind = 4 ) r4_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 4 ) x + + r1 = r4_uniform_01 ( seed ) + r2 = r4_uniform_01 ( seed ) + x = sqrt ( - 2.0E+00 * log ( r1 ) ) * cos ( 2.0E+00 * r4_pi * r2 ) + + i4_normal_ab = nint ( a + b * x ) + + return +end +function i8_normal_ab ( a, b, seed ) + +!*****************************************************************************80 +! +!! I8_NORMAL_AB returns a scaled pseudonormal I8. +! +! Discussion: +! +! The normal probability distribution function (PDF) is sampled, +! with mean A and standard deviation B. +! +! The result is then rounded to the nearest integer. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, real ( kind = 8 ) A, the mean of the PDF. +! +! Input, real ( kind = 8 ) B, the standard deviation of the PDF. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the +! random number generator. +! +! Output, integer ( kind = 8 ) I8_NORMAL_AB, a sample of the normal PDF. +! + implicit none + + real ( kind = 8 ) a + real ( kind = 8 ) b + integer ( kind = 8 ) i8_normal_ab + real ( kind = 8 ) r1 + real ( kind = 8 ) r2 + real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00 + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 8 ) x + + r1 = r8_uniform_01 ( seed ) + r2 = r8_uniform_01 ( seed ) + x = sqrt ( - 2.0D+00 * log ( r1 ) ) * cos ( 2.0D+00 * r8_pi * r2 ) + + i8_normal_ab = nint ( a + b * x ) + + return +end +function r4_normal_01 ( seed ) + +!*****************************************************************************80 +! +!! R4_NORMAL_01 returns a unit pseudonormal R4. +! +! Discussion: +! +! The standard normal probability distribution function (PDF) has +! mean 0 and standard deviation 1. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 4 ) R4_NORMAL_01, a sample of the standard +! normal PDF. +! + implicit none + + real ( kind = 4 ) r1 + real ( kind = 4 ) r2 + real ( kind = 4 ) r4_normal_01 + real ( kind = 4 ), parameter :: r4_pi = 3.141592653589793E+00 + real ( kind = 4 ) r4_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 4 ) x + + r1 = r4_uniform_01 ( seed ) + r2 = r4_uniform_01 ( seed ) + x = sqrt ( - 2.0E+00 * log ( r1 ) ) * cos ( 2.0E+00 * r4_pi * r2 ) + + r4_normal_01 = x + + return +end +function r4_normal_ab ( a, b, seed ) + +!*****************************************************************************80 +! +!! R4_NORMAL_AB returns a scaled pseudonormal R4. +! +! Discussion: +! +! The normal probability distribution function (PDF) is sampled, +! with mean A and standard deviation B. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, real ( kind = 4 ) A, the mean of the PDF. +! +! Input, real ( kind = 4 ) B, the standard deviation of the PDF. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 4 ) R4_NORMAL_AB, a sample of the normal PDF. +! + implicit none + + real ( kind = 4 ) a + real ( kind = 4 ) b + real ( kind = 4 ) r1 + real ( kind = 4 ) r2 + real ( kind = 4 ) r4_normal_ab + real ( kind = 4 ), parameter :: r4_pi = 3.141592653589793E+00 + real ( kind = 4 ) r4_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 4 ) x + + r1 = r4_uniform_01 ( seed ) + r2 = r4_uniform_01 ( seed ) + x = sqrt ( - 2.0E+00 * log ( r1 ) ) * cos ( 2.0E+00 * r4_pi * r2 ) + + r4_normal_ab = a + b * x + + return +end +function r4_uniform_01 ( seed ) + +!*****************************************************************************80 +! +!! R4_UNIFORM_01 returns a unit pseudorandom R4. +! +! Discussion: +! +! An R4 is a real ( kind = 4 ) value. +! +! This routine implements the recursion +! +! seed = 16807 * seed mod ( 2^31 - 1 ) +! r4_uniform_01 = seed / ( 2^31 - 1 ) +! +! The integer arithmetic never requires more than 32 bits, +! including a sign bit. +! +! If the initial seed is 12345, then the first three computations are +! +! Input Output R4_UNIFORM_01 +! SEED SEED +! +! 12345 207482415 0.096616 +! 207482415 1790989824 0.833995 +! 1790989824 2035175616 0.947702 +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! +! Reference: +! +! Paul Bratley, Bennett Fox, Linus Schrage, +! A Guide to Simulation, +! Second Edition, +! Springer, 1987, +! ISBN: 0387964673, +! LC: QA76.9.C65.B73. +! +! Bennett Fox, +! Algorithm 647: +! Implementation and Relative Efficiency of Quasirandom +! Sequence Generators, +! ACM Transactions on Mathematical Software, +! Volume 12, Number 4, December 1986, pages 362-376. +! +! Pierre L'Ecuyer, +! Random Number Generation, +! in Handbook of Simulation, +! edited by Jerry Banks, +! Wiley, 1998, +! ISBN: 0471134031, +! LC: T57.62.H37. +! +! Peter Lewis, Allen Goodman, James Miller, +! A Pseudo-Random Number Generator for the System/360, +! IBM Systems Journal, +! Volume 8, Number 2, 1969, pages 136-143. +! +! Parameters: +! +! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which +! should NOT be 0. On output, SEED has been updated. +! +! Output, real ( kind = 4 ) R4_UNIFORM_01, a new pseudorandom variate, +! strictly between 0 and 1. +! + implicit none + + integer ( kind = 4 ), parameter :: i4_huge = 2147483647 + integer ( kind = 4 ) k + integer ( kind = 4 ) seed + real ( kind = 4 ) r4_uniform_01 + + if ( seed == 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R4_UNIFORM_01 - Fatal error!' + write ( *, '(a)' ) ' Input value of SEED = 0.' + stop 1 + end if + + k = seed / 127773 + + seed = 16807 * ( seed - k * 127773 ) - k * 2836 + + if ( seed < 0 ) then + seed = seed + i4_huge + end if + + r4_uniform_01 = real ( seed, kind = 4 ) * 4.656612875E-10 + + return +end +subroutine r4vec_uniform_01 ( n, seed, r ) + +!*****************************************************************************80 +! +!! R4VEC_UNIFORM_01 returns a unit pseudorandom R4VEC. +! +! Discussion: +! +! An R4VEC is an array of real ( kind = 4 ) values. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! +! Reference: +! +! Paul Bratley, Bennett Fox, Linus Schrage, +! A Guide to Simulation, +! Second Edition, +! Springer, 1987, +! ISBN: 0387964673, +! LC: QA76.9.C65.B73. +! +! Bennett Fox, +! Algorithm 647: +! Implementation and Relative Efficiency of Quasirandom +! Sequence Generators, +! ACM Transactions on Mathematical Software, +! Volume 12, Number 4, December 1986, pages 362-376. +! +! Pierre L'Ecuyer, +! Random Number Generation, +! in Handbook of Simulation, +! edited by Jerry Banks, +! Wiley, 1998, +! ISBN: 0471134031, +! LC: T57.62.H37. +! +! Peter Lewis, Allen Goodman, James Miller, +! A Pseudo-Random Number Generator for the System/360, +! IBM Systems Journal, +! Volume 8, Number 2, 1969, pages 136-143. +! +! Parameters: +! +! Input, integer ( kind = 4 ) N, the number of entries in the vector. +! +! Input/output, integer ( kind = 4 ) SEED, the "seed" value, +! which should NOT be 0. +! On output, SEED has been updated. +! +! Output, real ( kind = 4 ) R(N), the vector of pseudorandom values. +! + implicit none + + integer ( kind = 4 ) n + + integer ( kind = 4 ) i + integer ( kind = 4 ), parameter :: i4_huge = 2147483647 + integer ( kind = 4 ) k + integer ( kind = 4 ) seed + real ( kind = 4 ) r(n) + + if ( seed == 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R4VEC_UNIFORM_01 - Fatal error!' + write ( *, '(a)' ) ' Input value of SEED = 0.' + stop 1 + end if + + do i = 1, n + + k = seed / 127773 + + seed = 16807 * ( seed - k * 127773 ) - k * 2836 + + if ( seed < 0 ) then + seed = seed + i4_huge + end if + + r(i) = real ( seed, kind = 4 ) * 4.656612875E-10 + + end do + + return +end +subroutine r4vec_normal_ab ( n, a, b, seed, x ) + +!*****************************************************************************80 +! +!! R4VEC_NORMAL_AB returns a scaled pseudonormal R4VEC. +! +! Discussion: +! +! The standard normal probability distribution function (PDF) has +! mean 0 and standard deviation 1. +! +! An R4VEC is a vector of R4's. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) N, the number of values desired. +! +! Input, real ( kind = 4 ) A, B, the mean and standard deviation. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 4 ) X(N), a sample of the standard normal PDF. +! +! Local parameters: +! +! Local, real ( kind = 4 ) R(N+1), is used to store some uniform +! random values. Its dimension is N+1, but really it is only needed +! to be the smallest even number greater than or equal to N. +! +! Local, integer ( kind = 4 ) X_LO_INDEX, X_HI_INDEX, records the range +! of entries of X that we need to compute. +! + implicit none + + integer ( kind = 4 ) n + + real ( kind = 4 ) a + real ( kind = 4 ) b + integer ( kind = 4 ) m + real ( kind = 4 ) r(n+1) + real ( kind = 4 ), parameter :: r4_pi = 3.141592653589793E+00 + real ( kind = 4 ) r4_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 4 ) x(n) + integer ( kind = 4 ) x_hi_index + integer ( kind = 4 ) x_lo_index +! +! Record the range of X we need to fill in. +! + x_lo_index = 1 + x_hi_index = n +! +! If we need just one new value, do that here to avoid null arrays. +! + if ( x_hi_index - x_lo_index + 1 == 1 ) then + + r(1) = r4_uniform_01 ( seed ) + + if ( r(1) == 0.0E+00 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R4VEC_NORMAL_AB - Fatal error!' + write ( *, '(a)' ) ' R4_UNIFORM_01 returned a value of 0.' + stop 1 + end if + + r(2) = r4_uniform_01 ( seed ) + + x(x_hi_index) = & + sqrt ( - 2.0E+00 * log ( r(1) ) ) * cos ( 2.0E+00 * r4_pi * r(2) ) +! +! If we require an even number of values, that's easy. +! + else if ( mod ( x_hi_index - x_lo_index, 2 ) == 1 ) then + + m = ( x_hi_index - x_lo_index + 1 ) / 2 + + call r4vec_uniform_01 ( 2*m, seed, r ) + + x(x_lo_index:x_hi_index-1:2) = & + sqrt ( - 2.0E+00 * log ( r(1:2*m-1:2) ) ) & + * cos ( 2.0E+00 * r4_pi * r(2:2*m:2) ) + + x(x_lo_index+1:x_hi_index:2) = & + sqrt ( - 2.0E+00 * log ( r(1:2*m-1:2) ) ) & + * sin ( 2.0E+00 * r4_pi * r(2:2*m:2) ) +! +! If we require an odd number of values, we generate an even number, +! and handle the last pair specially, storing one in X(N), and +! saving the other for later. +! + else + + x_hi_index = x_hi_index - 1 + + m = ( x_hi_index - x_lo_index + 1 ) / 2 + 1 + + call r4vec_uniform_01 ( 2*m, seed, r ) + + x(x_lo_index:x_hi_index-1:2) = & + sqrt ( - 2.0E+00 * log ( r(1:2*m-3:2) ) ) & + * cos ( 2.0E+00 * r4_pi * r(2:2*m-2:2) ) + + x(x_lo_index+1:x_hi_index:2) = & + sqrt ( - 2.0E+00 * log ( r(1:2*m-3:2) ) ) & + * sin ( 2.0E+00 * r4_pi * r(2:2*m-2:2) ) + + x(n) = sqrt ( - 2.0E+00 * log ( r(2*m-1) ) ) & + * cos ( 2.0E+00 * r4_pi * r(2*m) ) + + end if + + x(1:n) = a + b * x(1:n) + + return +end +function r8_normal_01 ( seed ) + +!*****************************************************************************80 +! +!! R8_NORMAL_01 returns a unit pseudonormal R8. +! +! Discussion: +! +! The standard normal probability distribution function (PDF) has +! mean 0 and standard deviation 1. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 8 ) R8_NORMAL_01, a normally distributed +! random value. +! + implicit none + + real ( kind = 8 ) r1 + real ( kind = 8 ) r2 + real ( kind = 8 ) r8_normal_01 + real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00 + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 8 ) x + + r1 = r8_uniform_01 ( seed ) + r2 = r8_uniform_01 ( seed ) + x = sqrt ( - 2.0D+00 * log ( r1 ) ) * cos ( 2.0D+00 * r8_pi * r2 ) + + r8_normal_01 = x + + return +end +function r8_normal_ab ( a, b, seed ) + +!*****************************************************************************80 +! +!! R8_NORMAL_AB returns a scaled pseudonormal R8. +! +! Discussion: +! +! The normal probability distribution function (PDF) is sampled, +! with mean A and standard deviation B. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, real ( kind = 8 ) A, the mean of the PDF. +! +! Input, real ( kind = 8 ) B, the standard deviation of the PDF. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 8 ) R8_NORMAL_AB, a sample of the normal PDF. +! + implicit none + + real ( kind = 8 ) a + real ( kind = 8 ) b + real ( kind = 8 ) r1 + real ( kind = 8 ) r2 + real ( kind = 8 ) r8_normal_ab + real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00 + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 8 ) x + + r1 = r8_uniform_01 ( seed ) + r2 = r8_uniform_01 ( seed ) + x = sqrt ( - 2.0D+00 * log ( r1 ) ) * cos ( 2.0D+00 * r8_pi * r2 ) + + r8_normal_ab = a + b * x + + return +end +function r8_uniform_01 ( seed ) + +!*****************************************************************************80 +! +!! R8_UNIFORM_01 returns a unit pseudorandom R8. +! +! Discussion: +! +! This routine implements the recursion +! +! seed = 16807 * seed mod ( 2^31 - 1 ) +! r8_uniform_01 = seed / ( 2^31 - 1 ) +! +! The integer arithmetic never requires more than 32 bits, +! including a sign bit. +! +! If the initial seed is 12345, then the first three computations are +! +! Input Output R8_UNIFORM_01 +! SEED SEED +! +! 12345 207482415 0.096616 +! 207482415 1790989824 0.833995 +! 1790989824 2035175616 0.947702 +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! +! Reference: +! +! Paul Bratley, Bennett Fox, Linus Schrage, +! A Guide to Simulation, +! Second Edition, +! Springer, 1987, +! ISBN: 0387964673, +! LC: QA76.9.C65.B73. +! +! Bennett Fox, +! Algorithm 647: +! Implementation and Relative Efficiency of Quasirandom +! Sequence Generators, +! ACM Transactions on Mathematical Software, +! Volume 12, Number 4, December 1986, pages 362-376. +! +! Pierre L'Ecuyer, +! Random Number Generation, +! in Handbook of Simulation, +! edited by Jerry Banks, +! Wiley, 1998, +! ISBN: 0471134031, +! LC: T57.62.H37. +! +! Peter Lewis, Allen Goodman, James Miller, +! A Pseudo-Random Number Generator for the System/360, +! IBM Systems Journal, +! Volume 8, 1969, pages 136-143. +! +! Parameters: +! +! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which +! should NOT be 0. +! On output, SEED has been updated. +! +! Output, real ( kind = 8 ) R8_UNIFORM_01, a new pseudorandom variate, +! strictly between 0 and 1. +! + implicit none + + integer ( kind = 8 ) k + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + + k = seed / 127773 + + seed = 16807 * ( seed - k * 127773 ) - k * 2836 + + if ( seed < 0 ) then + seed = seed + 2147483647 + end if +! +! Although SEED can be represented exactly as a 32 bit integer, +! it generally cannot be represented exactly as a 32 bit real number! +! + r8_uniform_01 = real ( seed, kind = 8 ) * 4.656612875D-10 + + return +end +subroutine r8mat_normal_01 ( m, n, seed, r ) + +!*****************************************************************************80 +! +!! R8MAT_NORMAL_01 returns a unit pseudonormal R8MAT. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 12 November 2010 +! +! Author: +! +! John Burkardt +! +! Reference: +! +! Paul Bratley, Bennett Fox, Linus Schrage, +! A Guide to Simulation, +! Second Edition, +! Springer, 1987, +! ISBN: 0387964673, +! LC: QA76.9.C65.B73. +! +! Bennett Fox, +! Algorithm 647: +! Implementation and Relative Efficiency of Quasirandom +! Sequence Generators, +! ACM Transactions on Mathematical Software, +! Volume 12, Number 4, December 1986, pages 362-376. +! +! Pierre L'Ecuyer, +! Random Number Generation, +! in Handbook of Simulation, +! edited by Jerry Banks, +! Wiley, 1998, +! ISBN: 0471134031, +! LC: T57.62.H37. +! +! Peter Lewis, Allen Goodman, James Miller, +! A Pseudo-Random Number Generator for the System/360, +! IBM Systems Journal, +! Volume 8, 1969, pages 136-143. +! +! Parameters: +! +! Input, integer ( kind = 4 ) M, N, the number of rows and columns +! in the array. +! +! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which +! should NOT be 0. On output, SEED has been updated. +! +! Output, real ( kind = 8 ) R(M,N), the array of pseudonormal values. +! + implicit none + + integer ( kind = 4 ) m + integer ( kind = 4 ) n + + integer ( kind = 4 ) seed + real ( kind = 8 ) r(m,n) + + call r8vec_normal_01 ( m * n, seed, r ) + + return +end +subroutine r8mat_normal_ab ( m, n, a, b, seed, r ) + +!*****************************************************************************80 +! +!! R8MAT_NORMAL_AB returns a scaled pseudonormal R8MAT. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Reference: +! +! Paul Bratley, Bennett Fox, Linus Schrage, +! A Guide to Simulation, +! Second Edition, +! Springer, 1987, +! ISBN: 0387964673, +! LC: QA76.9.C65.B73. +! +! Bennett Fox, +! Algorithm 647: +! Implementation and Relative Efficiency of Quasirandom +! Sequence Generators, +! ACM Transactions on Mathematical Software, +! Volume 12, Number 4, December 1986, pages 362-376. +! +! Pierre L'Ecuyer, +! Random Number Generation, +! in Handbook of Simulation, +! edited by Jerry Banks, +! Wiley, 1998, +! ISBN: 0471134031, +! LC: T57.62.H37. +! +! Peter Lewis, Allen Goodman, James Miller, +! A Pseudo-Random Number Generator for the System/360, +! IBM Systems Journal, +! Volume 8, 1969, pages 136-143. +! +! Parameters: +! +! Input, integer ( kind = 4 ) M, N, the number of rows and columns +! in the array. +! +! Input, real ( kind = 8 ) A, B, the mean and standard deviation. +! +! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which +! should NOT be 0. On output, SEED has been updated. +! +! Output, real ( kind = 8 ) R(M,N), the array of pseudonormal values. +! + implicit none + + integer ( kind = 4 ) m + integer ( kind = 4 ) n + + real ( kind = 8 ) a + real ( kind = 8 ) b + integer ( kind = 4 ) seed + real ( kind = 8 ) r(m,n) + + call r8vec_normal_ab ( m * n, a, b, seed, r ) + + return +end +subroutine r8mat_print ( m, n, a, title ) + +!*****************************************************************************80 +! +!! R8MAT_PRINT prints an R8MAT. +! +! Discussion: +! +! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 12 September 2004 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) M, the number of rows in A. +! +! Input, integer ( kind = 4 ) N, the number of columns in A. +! +! Input, real ( kind = 8 ) A(M,N), the matrix. +! +! Input, character ( len = * ) TITLE, a title. +! + implicit none + + integer ( kind = 4 ) m + integer ( kind = 4 ) n + + real ( kind = 8 ) a(m,n) + character ( len = * ) title + + call r8mat_print_some ( m, n, a, 1, 1, m, n, title ) + + return +end +subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) + +!*****************************************************************************80 +! +!! R8MAT_PRINT_SOME prints some of an R8MAT. +! +! Discussion: +! +! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 10 September 2009 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) M, N, the number of rows and columns. +! +! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. +! +! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. +! +! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. +! +! Input, character ( len = * ) TITLE, a title. +! + implicit none + + integer ( kind = 4 ), parameter :: incx = 5 + integer ( kind = 4 ) m + integer ( kind = 4 ) n + + real ( kind = 8 ) a(m,n) + character ( len = 14 ) ctemp(incx) + integer ( kind = 4 ) i + integer ( kind = 4 ) i2hi + integer ( kind = 4 ) i2lo + integer ( kind = 4 ) ihi + integer ( kind = 4 ) ilo + integer ( kind = 4 ) inc + integer ( kind = 4 ) j + integer ( kind = 4 ) j2 + integer ( kind = 4 ) j2hi + integer ( kind = 4 ) j2lo + integer ( kind = 4 ) jhi + integer ( kind = 4 ) jlo + character ( len = * ) title + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) trim ( title ) + + if ( m <= 0 .or. n <= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) ' (None)' + return + end if + + do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx + + j2hi = j2lo + incx - 1 + j2hi = min ( j2hi, n ) + j2hi = min ( j2hi, jhi ) + + inc = j2hi + 1 - j2lo + + write ( *, '(a)' ) ' ' + + do j = j2lo, j2hi + j2 = j + 1 - j2lo + write ( ctemp(j2), '(i8,6x)' ) j + end do + + write ( *, '('' Col '',5a14)' ) ctemp(1:inc) + write ( *, '(a)' ) ' Row' + write ( *, '(a)' ) ' ' + + i2lo = max ( ilo, 1 ) + i2hi = min ( ihi, m ) + + do i = i2lo, i2hi + + do j2 = 1, inc + + j = j2lo - 1 + j2 + + if ( a(i,j) == real ( int ( a(i,j) ), kind = 8 ) ) then + write ( ctemp(j2), '(f8.0,6x)' ) a(i,j) + else + write ( ctemp(j2), '(g14.6)' ) a(i,j) + end if + + end do + + write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc ) + + end do + + end do + + return +end +subroutine r8vec_normal_01 ( n, seed, x ) + +!*****************************************************************************80 +! +!! R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC. +! +! Discussion: +! +! An R8VEC is an array of double precision real values. +! +! The standard normal probability distribution function (PDF) has +! mean 0 and standard deviation 1. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 18 May 2014 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) N, the number of values desired. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 8 ) X(N), a sample of the standard normal PDF. +! +! Local parameters: +! +! Local, real ( kind = 8 ) R(N+1), is used to store some uniform +! random values. Its dimension is N+1, but really it is only needed +! to be the smallest even number greater than or equal to N. +! +! Local, integer ( kind = 4 ) X_LO_INDEX, X_HI_INDEX, records the range +! of entries of X that we need to compute +! + implicit none + + integer ( kind = 4 ) n + + integer ( kind = 4 ) m + real ( kind = 8 ) r(n+1) + real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00 + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 8 ) x(n) + integer ( kind = 4 ) x_hi_index + integer ( kind = 4 ) x_lo_index +! +! Record the range of X we need to fill in. +! + x_lo_index = 1 + x_hi_index = n +! +! If we need just one new value, do that here to avoid null arrays. +! + if ( x_hi_index - x_lo_index + 1 == 1 ) then + + r(1) = r8_uniform_01 ( seed ) + + if ( r(1) == 0.0D+00 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8VEC_NORMAL_01 - Fatal error!' + write ( *, '(a)' ) ' R8_UNIFORM_01 returned a value of 0.' + stop 1 + end if + + r(2) = r8_uniform_01 ( seed ) + + x(x_hi_index) = & + sqrt ( - 2.0D+00 * log ( r(1) ) ) * cos ( 2.0D+00 * r8_pi * r(2) ) +! +! If we require an even number of values, that's easy. +! + else if ( mod ( x_hi_index - x_lo_index, 2 ) == 1 ) then + + m = ( x_hi_index - x_lo_index + 1 ) / 2 + + call r8vec_uniform_01 ( 2*m, seed, r ) + + x(x_lo_index:x_hi_index-1:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-1:2) ) ) & + * cos ( 2.0D+00 * r8_pi * r(2:2*m:2) ) + + x(x_lo_index+1:x_hi_index:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-1:2) ) ) & + * sin ( 2.0D+00 * r8_pi * r(2:2*m:2) ) +! +! If we require an odd number of values, we generate an even number, +! and handle the last pair specially, storing one in X(N), and +! saving the other for later. +! + else + + x_hi_index = x_hi_index - 1 + + m = ( x_hi_index - x_lo_index + 1 ) / 2 + 1 + + call r8vec_uniform_01 ( 2*m, seed, r ) + + x(x_lo_index:x_hi_index-1:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-3:2) ) ) & + * cos ( 2.0D+00 * r8_pi * r(2:2*m-2:2) ) + + x(x_lo_index+1:x_hi_index:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-3:2) ) ) & + * sin ( 2.0D+00 * r8_pi * r(2:2*m-2:2) ) + + x(n) = sqrt ( - 2.0D+00 * log ( r(2*m-1) ) ) & + * cos ( 2.0D+00 * r8_pi * r(2*m) ) + + end if + + return +end +subroutine r8vec_normal_ab ( n, a, b, seed, x ) + +!*****************************************************************************80 +! +!! R8VEC_NORMAL_AB returns a scaled pseudonormal R8VEC. +! +! Discussion: +! +! The standard normal probability distribution function (PDF) has +! mean 0 and standard deviation 1. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 06 August 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) N, the number of values desired. +! +! Input, real ( kind = 8 ) A, B, the mean and standard deviation. +! +! Input/output, integer ( kind = 4 ) SEED, a seed for the random +! number generator. +! +! Output, real ( kind = 8 ) X(N), a sample of the standard normal PDF. +! +! Local parameters: +! +! Local, real ( kind = 8 ) R(N+1), is used to store some uniform +! random values. Its dimension is N+1, but really it is only needed +! to be the smallest even number greater than or equal to N. +! +! Local, integer ( kind = 4 ) X_LO_INDEX, X_HI_INDEX, records the range +! of entries of X that we need to compute. +! + implicit none + + integer ( kind = 4 ) n + + real ( kind = 8 ) a + real ( kind = 8 ) b + integer ( kind = 4 ) m + real ( kind = 8 ) r(n+1) + real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00 + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + real ( kind = 8 ) x(n) + integer ( kind = 4 ) x_hi_index + integer ( kind = 4 ) x_lo_index +! +! Record the range of X we need to fill in. +! + x_lo_index = 1 + x_hi_index = n +! +! If we need just one new value, do that here to avoid null arrays. +! + if ( x_hi_index - x_lo_index + 1 == 1 ) then + + r(1) = r8_uniform_01 ( seed ) + + if ( r(1) == 0.0D+00 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8VEC_NORMAL_AB - Fatal error!' + write ( *, '(a)' ) ' R8_UNIFORM_01 returned a value of 0.' + stop 1 + end if + + r(2) = r8_uniform_01 ( seed ) + + x(x_hi_index) = & + sqrt ( - 2.0D+00 * log ( r(1) ) ) * cos ( 2.0D+00 * r8_pi * r(2) ) +! +! If we require an even number of values, that's easy. +! + else if ( mod ( x_hi_index - x_lo_index, 2 ) == 1 ) then + + m = ( x_hi_index - x_lo_index + 1 ) / 2 + + call r8vec_uniform_01 ( 2*m, seed, r ) + + x(x_lo_index:x_hi_index-1:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-1:2) ) ) & + * cos ( 2.0D+00 * r8_pi * r(2:2*m:2) ) + + x(x_lo_index+1:x_hi_index:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-1:2) ) ) & + * sin ( 2.0D+00 * r8_pi * r(2:2*m:2) ) +! +! If we require an odd number of values, we generate an even number, +! and handle the last pair specially, storing one in X(N), and +! saving the other for later. +! + else + + x_hi_index = x_hi_index - 1 + + m = ( x_hi_index - x_lo_index + 1 ) / 2 + 1 + + call r8vec_uniform_01 ( 2*m, seed, r ) + + x(x_lo_index:x_hi_index-1:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-3:2) ) ) & + * cos ( 2.0D+00 * r8_pi * r(2:2*m-2:2) ) + + x(x_lo_index+1:x_hi_index:2) = & + sqrt ( - 2.0D+00 * log ( r(1:2*m-3:2) ) ) & + * sin ( 2.0D+00 * r8_pi * r(2:2*m-2:2) ) + + x(n) = sqrt ( - 2.0D+00 * log ( r(2*m-1) ) ) & + * cos ( 2.0D+00 * r8_pi * r(2*m) ) + + end if + + x(1:n) = a + b * x(1:n) + + return +end +subroutine r8vec_print ( n, a, title ) + +!*****************************************************************************80 +! +!! R8VEC_PRINT prints an R8VEC. +! +! Discussion: +! +! An R8VEC is a vector of R8's. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 22 August 2000 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, integer ( kind = 4 ) N, the number of components of the vector. +! +! Input, real ( kind = 8 ) A(N), the vector to be printed. +! +! Input, character ( len = * ) TITLE, a title. +! + implicit none + + integer ( kind = 4 ) n + + real ( kind = 8 ) a(n) + integer ( kind = 4 ) i + character ( len = * ) title + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) trim ( title ) + write ( *, '(a)' ) ' ' + + do i = 1, n + write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i) + end do + + return +end +subroutine r8vec_uniform_01 ( n, seed, r ) + +!*****************************************************************************80 +! +!! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! +! Reference: +! +! Paul Bratley, Bennett Fox, Linus Schrage, +! A Guide to Simulation, +! Second Edition, +! Springer, 1987, +! ISBN: 0387964673, +! LC: QA76.9.C65.B73. +! +! Bennett Fox, +! Algorithm 647: +! Implementation and Relative Efficiency of Quasirandom +! Sequence Generators, +! ACM Transactions on Mathematical Software, +! Volume 12, Number 4, December 1986, pages 362-376. +! +! Pierre L'Ecuyer, +! Random Number Generation, +! in Handbook of Simulation, +! edited by Jerry Banks, +! Wiley, 1998, +! ISBN: 0471134031, +! LC: T57.62.H37. +! +! Peter Lewis, Allen Goodman, James Miller, +! A Pseudo-Random Number Generator for the System/360, +! IBM Systems Journal, +! Volume 8, 1969, pages 136-143. +! +! Parameters: +! +! Input, integer ( kind = 4 ) N, the number of entries in the vector. +! +! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which +! should NOT be 0. On output, SEED has been updated. +! +! Output, real ( kind = 8 ) R(N), the vector of pseudorandom values. +! + implicit none + + integer ( kind = 4 ) n + + integer ( kind = 4 ) i + integer ( kind = 4 ) k + integer ( kind = 4 ) seed + real ( kind = 8 ) r(n) + + do i = 1, n + + k = seed / 127773 + + seed = 16807 * ( seed - k * 127773 ) - k * 2836 + + if ( seed < 0 ) then + seed = seed + 2147483647 + end if + + r(i) = real ( seed, kind = 8 ) * 4.656612875D-10 + + end do + + return +end +subroutine timestamp ( ) + +!*****************************************************************************80 +! +!! TIMESTAMP prints the current YMDHMS date as a time stamp. +! +! Example: +! +! 31 May 2001 9:45:54.872 AM +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 18 May 2013 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! None +! + implicit none + + character ( len = 8 ) ampm + integer ( kind = 4 ) d + integer ( kind = 4 ) h + integer ( kind = 4 ) m + integer ( kind = 4 ) mm + character ( len = 9 ), parameter, dimension(12) :: month = (/ & + 'January ', 'February ', 'March ', 'April ', & + 'May ', 'June ', 'July ', 'August ', & + 'September', 'October ', 'November ', 'December ' /) + integer ( kind = 4 ) n + integer ( kind = 4 ) s + integer ( kind = 4 ) values(8) + integer ( kind = 4 ) y + + call date_and_time ( values = values ) + + y = values(1) + m = values(2) + d = values(3) + h = values(5) + n = values(6) + s = values(7) + mm = values(8) + + if ( h < 12 ) then + ampm = 'AM' + else if ( h == 12 ) then + if ( n == 0 .and. s == 0 ) then + ampm = 'Noon' + else + ampm = 'PM' + end if + else + h = h - 12 + if ( h < 12 ) then + ampm = 'PM' + else if ( h == 12 ) then + if ( n == 0 .and. s == 0 ) then + ampm = 'Midnight' + else + ampm = 'AM' + end if + end if + end if + + write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & + d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) + + return +end diff --git a/normal_prb.f90 b/normal_prb.f90 new file mode 100644 index 0000000..0301b8d --- /dev/null +++ b/normal_prb.f90 @@ -0,0 +1,673 @@ +program main + +!*****************************************************************************80 +! +!! MAIN is the main program for NORMAL_PRB. +! +! Discussion: +! +! NORMAL_PRB tests the NORMAL library. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 03 March 2015 +! +! Author: +! +! John Burkardt +! + implicit none + + call timestamp ( ) + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'NORMAL_PRB' + write ( *, '(a)' ) ' FORTRAN90 version;' + write ( *, '(a)' ) ' Test the NORMAL library.' + + call c4_normal_01_test ( ) + call c8_normal_01_test ( ) + call i4_normal_ab_test ( ) + call r4_normal_01_test ( ) + call r4_normal_ab_test ( ) + call r4_uniform_01_test ( ) + call r8_normal_01_test ( ) + call r8_normal_ab_test ( ) + call r8_uniform_01_test ( ) + call r8mat_normal_01_test ( ) + call r8mat_normal_ab_test ( ) + call r8vec_normal_01_test ( ) + call r8vec_normal_ab_test ( ) + call r8vec_uniform_01_test ( ) +! +! Terminate. +! + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'NORMAL_PRB' + write ( *, '(a)' ) ' Normal end of execution.' + write ( *, '(a)' ) ' ' + call timestamp ( ) + + stop +end +subroutine c4_normal_01_test ( ) + +!*****************************************************************************80 +! +!! C4_NORMAL_01_TEST tests C4_NORMAL_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + complex ( kind = 4 ) c4_normal_01 + integer ( kind = 4 ) i + complex ( kind = 4 ) r + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'C4_NORMAL_01_TEST' + write ( *, '(a)' ) ' C4_NORMAL_01 computes pseudorandom complex values ' + write ( *, '(a)' ) ' with a Normal 01 circular distribution.' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = c4_normal_01 ( seed ) + write ( *, '(2x,i8,2x,f14.8,2x,f14.8)' ) i, r + end do + + return +end +subroutine c8_normal_01_test ( ) + +!*****************************************************************************80 +! +!! C8_NORMAL_01_TEST tests C8_NORMAL_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + complex ( kind = 8 ) c8_normal_01 + integer ( kind = 4 ) i + complex ( kind = 8 ) r + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'C8_NORMAL_01_TEST' + write ( *, '(a)' ) ' C8_NORMAL_01 computes pseudorandom double precision' + write ( *, '(a)' ) ' complex values with Normal 01 circular distribution.' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = c8_normal_01 ( seed ) + write ( *, '(2x,i8,2x,f14.8,2x,f14.8)' ) i, r + end do + + return +end +subroutine i4_normal_ab_test ( ) + +!*****************************************************************************80 +! +!! I4_NORMAL_AB_TEST tests I4_NORMAL_AB. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + integer ( kind = 4 ) i4_normal_ab + real ( kind = 4 ) mu + integer ( kind = 4 ) r + integer ( kind = 4 ) seed + real ( kind = 4 ) sigma + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'I4_NORMAL_AB_TEST' + write ( *, '(a)' ) ' I4_NORMAL_AB computes integer pseudonormal values ' + write ( *, '(a)' ) ' with mean MU and standard deviation SIGMA.' + + mu = 70.0E+00 + sigma = 10.0E+00 + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,g14.6)' ) ' MU = ', mu + write ( *, '(a,g14.6)' ) ' SIGMA = ', sigma + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = i4_normal_ab ( mu, sigma, seed ) + write ( *, '(2x,i8,2x,i8)' ) i, r + end do + + return +end +subroutine r4_normal_01_test ( ) + +!*****************************************************************************80 +! +!! R4_NORMAL_01_TEST tests R4_NORMAL_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + real ( kind = 4 ) r + real ( kind = 4 ) r4_normal_01 + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R4_NORMAL_01_TEST' + write ( *, '(a)' ) ' R4_NORMAL_01 computes normal pseudorandom values ' + write ( *, '(a)' ) ' with mean 0.0 and standard deviation 1.0.' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = r4_normal_01 ( seed ) + write ( *, '(2x,i8,2x,f14.8)' ) i, r + end do + + return +end +subroutine r4_normal_ab_test ( ) + +!*****************************************************************************80 +! +!! R4_NORMAL_AB_TEST tests R4_NORMAL_AB. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 03 March 2015 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + real ( kind = 4 ) mu + real ( kind = 4 ) r + real ( kind = 4 ) r4_normal_ab + integer ( kind = 4 ) seed + real ( kind = 4 ) sigma + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R4_NORMAL_AB_TEST' + write ( *, '(a)' ) ' R4_NORMAL_AB computes real pseudonormal values ' + write ( *, '(a)' ) ' with mean MU and standard deviation SIGMA.' + + mu = 10.0E+00 + sigma = 2.0E+00 + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,g14.6)' ) ' MU = ', mu + write ( *, '(a,g14.6)' ) ' SIGMA = ', sigma + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = r4_normal_ab ( mu, sigma, seed ) + write ( *, '(2x,i8,2x,f14.8)' ) i, r + end do + + return +end +subroutine r4_uniform_01_test ( ) + +!*****************************************************************************80 +! +!! R4_UNIFORM_01_TEST tests R4_UNIFORM_01 +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 03 March 2015 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + integer ( kind = 4 ) seed + real ( kind = 4 ) r + real ( kind = 4 ) r4_uniform_01 + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R4_UNIFORM_01_TEST' + write ( *, '(a)' ) ' R4_UNIFORM_01 samples a uniform random' + write ( *, '(a)' ) ' distribution in [0,1].' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) '' + + do i = 1, 10 + r = r4_uniform_01 ( seed ) + write ( *, '(2x,i2,2x,g14.6)' ) i, r + end do + + return +end +subroutine r8_normal_01_test ( ) + +!*****************************************************************************80 +! +!! R8_NORMAL_01_TEST tests R8_NORMAL_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + real ( kind = 8 ) r + real ( kind = 8 ) r8_normal_01 + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8_NORMAL_01_TEST' + write ( *, '(a)' ) ' R8_NORMAL_01 computes pseudonormal values ' + write ( *, '(a)' ) ' with mean 0.0 and standard deviation 1.0.' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = r8_normal_01 ( seed ) + write ( *, '(2x,i8,2x,g14.6)' ) i, r + end do + + return +end +subroutine r8_normal_ab_test ( ) + +!*****************************************************************************80 +! +!! R8_NORMAL_AB_TEST tests R8_NORMAL_AB. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + real ( kind = 8 ) mu + real ( kind = 8 ) r + real ( kind = 8 ) r8_normal_ab + integer ( kind = 4 ) seed + real ( kind = 8 ) sigma + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8_NORMAL_AB_TEST' + write ( *, '(a)' ) ' R8_NORMAL_AB computes pseudonormal values ' + write ( *, '(a)' ) ' with mean MU and standard deviation SIGMA.' + + mu = 10.0D+00 + sigma = 2.0D+00 + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,g14.6)' ) ' MU = ', mu + write ( *, '(a,g14.6)' ) ' SIGMA = ', sigma + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) ' ' + + do i = 1, 10 + r = r8_normal_ab ( mu, sigma, seed ) + write ( *, '(2x,i8,2x,f14.8)' ) i, r + end do + + return +end +subroutine r8_uniform_01_test ( ) + +!*****************************************************************************80 +! +!! R8_UNIFORM_01_TEST tests R8_UNIFORM_01 +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 05 December 2006 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ) i + real ( kind = 8 ) r + real ( kind = 8 ) r8_uniform_01 + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8_UNIFORM_01_TEST' + write ( *, '(a)' ) ' R8_UNIFORM_01 samples a uniform random' + write ( *, '(a)' ) ' distribution in [0,1].' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + write ( *, '(a)' ) '' + + do i = 1, 10 + r = r8_uniform_01 ( seed ) + write ( *, '(2x,i2,2x,g14.6)' ) i, r + end do + + return +end +subroutine r8mat_normal_01_test ( ) + +!*****************************************************************************80 +! +!! R8MAT_NORMAL_01_TEST tests R8MAT_NORMAL_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ), parameter :: m = 5 + integer ( kind = 4 ), parameter :: n = 4 + + real ( kind = 8 ) r(m,n) + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8MAT_NORMAL_01_TEST' + write ( *, '(a)' ) ' R8MAT_NORMAL_01 returns a matrix of Normal 01 values.' + + seed = 123456789 + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + + call r8mat_normal_01 ( m, n, seed, r ) + + call r8mat_print ( m, n, r, ' Matrix of Normal 01 values:' ) + + return +end +subroutine r8mat_normal_ab_test ( ) + +!*****************************************************************************80 +! +!! R8MAT_NORMAL_AB_TEST tests R8MAT_NORMAL_AB. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ), parameter :: m = 5 + integer ( kind = 4 ), parameter :: n = 4 + + real ( kind = 8 ) mu + real ( kind = 8 ) r(m,n) + integer ( kind = 4 ) seed + real ( kind = 8 ) sigma + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8MAT_NORMAL_AB_TEST' + write ( *, '(a)' ) ' R8MAT_NORMAL_AB returns a matrix of Normal AB values.' + + mu = 100.0D+00 + sigma = 5.0D+00 + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,g14.6)' ) ' MU = ', mu + write ( *, '(a,g14.6)' ) ' SIGMA = ', sigma + write ( *, '(a,i12)' ) ' SEED = ', seed + + call r8mat_normal_ab ( m, n, mu, sigma, seed, r ) + + call r8mat_print ( m, n, r, ' Matrix of Normal AB values:' ) + + return +end +subroutine r8vec_normal_01_test ( ) + +!*****************************************************************************80 +! +!! R8VEC_NORMAL_01_TEST tests R8VEC_NORMAL_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 31 May 2007 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ), parameter :: n = 10 + + real ( kind = 8 ) r(n) + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8VEC_NORMAL_01_TEST' + write ( *, '(a)' ) ' R8VEC_NORMAL_01 computes a vector of Normal 01 values.' + + seed = 123456789 + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + + call r8vec_normal_01 ( n, seed, r ) + + call r8vec_print ( n, r, ' Vector of Normal 01 values:' ) + + return +end +subroutine r8vec_normal_ab_test ( ) + +!*****************************************************************************80 +! +!! R8VEC_NORMAL_AB_TEST tests R8VEC_NORMAL_AB. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 03 March 2015 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ), parameter :: n = 10 + + real ( kind = 8 ) mu + real ( kind = 8 ) r(n) + integer ( kind = 4 ) seed + real ( kind = 8 ) sigma + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8VEC_NORMAL_AB_TEST' + write ( *, '(a)' ) ' R8VEC_NORMAL_AB computes a vector of Normal AB values.' + + mu = 15.0D+00 + sigma = 0.25D+00 + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,g14.6)' ) ' MU = ', mu + write ( *, '(a,g14.6)' ) ' SIGMA = ', sigma + write ( *, '(a,i12)' ) ' SEED = ', seed + + call r8vec_normal_ab ( n, mu, sigma, seed, r ) + + call r8vec_print ( n, r, ' Vector of Normal AB values:' ) + + return +end +subroutine r8vec_uniform_01_test ( ) + +!*****************************************************************************80 +! +!! R8VEC_UNIFORM_01_TEST tests R8VEC_UNIFORM_01. +! +! Licensing: +! +! This code is distributed under the GNU LGPL license. +! +! Modified: +! +! 25 June 2010 +! +! Author: +! +! John Burkardt +! + implicit none + + integer ( kind = 4 ), parameter :: n = 10 + + real ( kind = 8 ) r(n) + integer ( kind = 4 ) seed + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'R8VEC_UNIFORM_01_TEST' + write ( *, '(a)' ) ' R8VEC_UNIFORM_01 returns a random R8VEC ' + write ( *, '(a)' ) ' with entries in [0,1].' + + seed = 123456789 + + write ( *, '(a)' ) ' ' + write ( *, '(a,i12)' ) ' SEED = ', seed + + call r8vec_uniform_01 ( n, seed, r ) + + call r8vec_print ( n, r, ' Random R8VEC:' ) + + return +end diff --git a/particle_gen.py b/particle_gen.py new file mode 100644 index 0000000..cfab619 --- /dev/null +++ b/particle_gen.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Jan 27 22:47:50 2019 + +Este programa é para colocar as partículas em seus lugares iniciais + +@author: gabriel +""" + +import numpy as np +import matplotlib.pyplot as plt +from random import randint + +arquivo1 = 'grupo1.csv' +arquivo2 = 'grupo2.csv' + +N1 = 1600 +N2 = 6400 +# Posição das partículas +p1 = np.zeros([N1,2]) +p2 = np.zeros([N2,2]) + + +d1 = 160 +d2 = 40 + + +spac = 2**(1/6) + +#posição de referência +refe = np.array([97,100]) +cont = 0 + + + +for i in range(d2): + for j in range(d2): + p1[cont,:] = [i*spac+refe[0],j*spac+refe[1]] + cont = cont+1 + +refe = [30,50] +cont = 0 +for i in range(d1): + for j in range(d2): + p2[cont,:] = [i*spac+refe[0],j*spac+refe[1]] + cont = cont+1 + +plt.scatter(p1[:,0],p1[:,1], label='grupo1') +plt.scatter(p2[:,0],p2[:,1], label='grupo2') +plt.legend() +plt.show() + +with open(arquivo1,'w') as file: + for i in range(N1): + file.write('{},{}\n'.format(p1[i,0],p1[i,1])) + +with open(arquivo2,'w') as file: + for i in range(N2): + file.write('{},{}\n'.format(p2[i,0],p2[i,1])) + + + + +##################################### diff --git a/particle_gen2.py b/particle_gen2.py new file mode 100644 index 0000000..c0c01ee --- /dev/null +++ b/particle_gen2.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Jan 27 22:47:50 2019 + +Este programa é para colocar as partículas em seus lugares iniciais + +@author: gabriel +""" + +import numpy as np +import matplotlib.pyplot as plt +from random import randint + +arquivo1 = 'basin1.csv' +arquivo2 = 'basin2.csv' + +N1 = 395 +N2 = 17238 +# Posição das partículas +p1 = np.zeros([N1,2]) +p2 = np.zeros([N2,2]) + +d1 = 221 +d2 = 78 + +spac = 2**(1/6) + +#posição de referência +refe = np.array([125,150]) +cont = 0 + +# bola +p1[cont,:] = [refe[0],refe[1]] + +cont +=1 +r = spac +while N1 > cont: + n = int(np.pi*r/spac)*2 + for i in range(n): + p1[cont,:] = [r*np.cos(2*np.pi*i/n) + refe[0],r*np.sin(2*np.pi*i/n) + refe[1]] + cont += 1 + if cont == N1: + break + r += spac + + +# fundo + +refe = [.1,.1] +cont = 0 +for i in range(d1): + for j in range(d2): + p2[cont,:] = [i*spac+refe[0] ,j*spac+refe[1]] + cont = cont+1 + +plt.scatter(p1[:,0],p1[:,1], label='grupo1') +plt.scatter(p2[:,0],p2[:,1], label='grupo2') +plt.legend() +plt.show() + +with open(arquivo1,'w') as file: + for i in range(N1): + file.write('{},{}\n'.format(p1[i,0],p1[i,1])) + +with open(arquivo2,'w') as file: + for i in range(N2): + file.write('{},{}\n'.format(p2[i,0],p2[i,1])) + + + + +##################################### diff --git a/particle_gen3.py b/particle_gen3.py new file mode 100644 index 0000000..0e881e3 --- /dev/null +++ b/particle_gen3.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Jan 27 22:47:50 2019 + +Este programa é para colocar as partículas em seus lugares iniciais + +@author: gabriel +""" + +import numpy as np +import matplotlib.pyplot as plt +from random import randint + +arquivo1 = 'rt1.csv' +arquivo2 = 'rt2.csv' + + +d1 = 155 +d2 = 20 + +N1 = d1*d2 +N2 = d1*d2 +# Posição das partículas +p1 = np.zeros([N1,2]) +p2 = np.zeros([N2,2]) + + +spac = 2**(1/6)*.8 + +#posição de referência +refe = np.array([0.1,37.4]) +cont = 0 + + + +for i in range(d1): + for j in range(d2): + p1[cont,:] = [i*spac+refe[0],-j*spac+refe[1]] + cont = cont+1 + +refe = [0.1,0.1] +cont = 0 +for i in range(d1): + for j in range(d2): + p2[cont,:] = [i*spac+refe[0],j*spac+refe[1]] + cont = cont+1 + +plt.scatter(p1[:,0],p1[:,1], label='rt1') +plt.scatter(p2[:,0],p2[:,1], label='rt2') +plt.legend() +plt.show() + +with open(arquivo1,'w') as file: + for i in range(N1): + file.write('{},{}\n'.format(p1[i,0],p1[i,1])) + +with open(arquivo2,'w') as file: + for i in range(N2): + file.write('{},{}\n'.format(p2[i,0],p2[i,1])) + + + + +##################################### diff --git a/randnormal.f90 b/randnormal.f90 new file mode 100644 index 0000000..f339485 --- /dev/null +++ b/randnormal.f90 @@ -0,0 +1,35 @@ +module randnormal + contains + + function GaussDeviate() + ! BOX MULLER + real :: GaussDeviate + real :: a1, a2, s, r, b1 + logical, save :: iset = .false. + real, save :: b2 + + r = 3.0 + + if (.not. iset) then + do while (r >= 1) + ! gerar números aleatórios entre -1 e 1 + call random_seed() + call random_number(a1) + call random_number(a2) + a1 = 2*a1-1 + a2 = 2*a2-1 + r = a1**2 + a2**2 + end do + s = sqrt(-2*log(r)/r) !transformação polar + b1 = a1*s + b2 = a2*s + iset = .true. + GaussDeviate = b1 + else + iset = .false. + GaussDeviate = b2 + end if + + end function + +end module randnormal diff --git a/settings.ini b/settings.ini new file mode 100644 index 0000000..64bf1ff --- /dev/null +++ b/settings.ini @@ -0,0 +1,104 @@ +# PARAMETROS GLOBAIS +[global] + nimpre = 150 # quntidade de saidas + N = 6200 # número de partículas + Ntype = 2 # número de tipos de partícula + t_fim = 14 # tempo de execução + dt = 0.0005 # passo de tempo + dimX = 140 # Dimensões da região de cálculo + dimY = 37.5 # + mesh = 56 15 # malha(elementos) + rcut = 2.5 # *sigma + wall = 'eeee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste + Td = 2.8971e+24 #temperatura desejada constante (-1 = não aplicar) + temp_Td = 0 14 1000 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] + vd = 0.1 0.1 # velocidade distribuida com Maxwell Boltzmann + fric_term = 0 # fricção artifical + 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 -12.44 # Campo de força gravitacional que afeta todas as partículas + +# PARTICLE PROPERTIES +[par_0] + x = 'rt1.csv' + v = 0 0 + m = 2 + nome = 'g1' + sigma = 0.9412 # Lennard Jones + epsilon = 1 # Lennard Jones + quantidade = 3100 + x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay + +[par_1] + x = 'rt2.csv' + v = 0 0 + m = 1 + nome = 'g2' + sigma = 1 # Lennard Jones + epsilon = 1 # Lennard Jones + quantidade = 3100 + x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay + +# Sobre NMPT. Pode-se estimar dividindo o tamanho da fronteira do domínio do processo por 2.5^(1/6) ou pelo raio. E multiplicar por +# um fator. + +#[par_1] +# x = 'grupo2.csv' +# v = 0 0 +# m = 1 +# nome = 'g2' +# quantidade = 2 +# sigma = 1 # Leonard Jones +# epsilon = 5 # Leonard Jones +# +#[par_2] +# x = 'grupo7.csv' +# v = -10 10 +# m = 1 +# nome = 'g3' +# quantidade = 8 +# +#[par_3] +# x = 'grupo8.csv' +# v = 0 0 +# m = 1 +# nome = 'g4' +# quantidade = 8 +# +#[par_0] +# x = 'grupo0.csv' +# v = 30 5 +# m = 1 +# nome = 'g1' +# quantidade = 2 +# sigma = 1 # Leonard Jones +# epsilon = 5 # Leonard Jones + +#nimpre = 150 # quantidade de saidas +#N = 1600 # número de partículas +#Ntype = 1 # número de tipos de partícula +#t_fim = 1 # tempo de execução +#dt = 0.00005 # passo de tempo +#dimX = 250 # Dimensões da região de cálculo +#dimY = 250 # +#mesh = 1 1 # malha(elementos) +#sigma = 1 # Leonard Jones +#epsilon = 5 # Leonard Jones +#rcut = 2.5 # *sigma +#wall = 'eeee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste +#Td = -1 +#fric_term = 0 +# +#[par_0] +# x = 'grupo1.csv' +# v = 0 -10 +# m = 1 +# nome = 'g1' +# quantidade = 1600 +# +#[par_1] +# x = 'grupo2.csv' +# v = 0 0 +# m = 1 +# nome = 'g2' +# quantidade = 6400 + diff --git a/teste.f90 b/teste.f90 new file mode 100644 index 0000000..4d90842 --- /dev/null +++ b/teste.f90 @@ -0,0 +1,12 @@ +program teste + use mod1 + + call system('mkdir temp') + open(10,file='temp/test.txt',status="replace") + write(10,*) 'lalalala' + print*,'a' + + + + +end program teste diff --git a/teste2.f90 b/teste2.f90 new file mode 100644 index 0000000..5b71049 --- /dev/null +++ b/teste2.f90 @@ -0,0 +1,504 @@ +! A derived type for storing data. +module mod1 + integer, parameter:: dp=kind(0.d0) ! double precision + ! vetor de strings + type string + character(len=:), allocatable :: str + integer :: quanty + end type string +end module mod1 + +module data + use linkedlist + implicit none + + private + public :: data_t + public :: data_ptr + + ! Data is stored in data_t + type :: data_t + real :: x + end type data_t + + ! A trick to allow us to store pointers in the list + type :: data_ptr + type(data_t), pointer :: p + end type data_ptr + +end module data + +module mod3 + use linkedlist + type :: container + type(list_t), pointer :: list + end type container +end module mod3 + +module matprint + contains + subroutine printm(m,r) + use mod1 + real(dp), dimension(:,:), intent(in) :: m + integer,dimension(2) :: s + integer :: i,j + character(*), optional, intent(in) :: r + + s = shape(m) + if(present(r)) then + print*, r,' =' + else + print*, ' =' + end if + + do i = 1,s(1) + do j = 1,s(2) + write(*,'(EN14.3)',advance='no') m(i,j) + end do + write(*,*) + end do + write(*,*) ' ' + end subroutine printm +end module matprint + +module mod2 + contains + + subroutine testes(m,ptr,ptrn) + use linkedlist + use mod3 + type(container),intent(in), allocatable,dimension(:,:) :: m + integer,pointer :: ptr,ptrn + print*,'linha 72' + ptrn = transfer(list_get(m(2,1)%list), ptr) + print *, 'Head node data:', ptrn + end subroutine testes + + subroutine testchange(node2,m) + use linkedlist + use mod3 + type(container),intent(in), allocatable,dimension(:,:) :: m + type(list_t), pointer :: node2 + call list_change(node2,m(2,1)%list) + end subroutine testchange + + subroutine sub2(y) + integer, intent(inout) :: y + + y = y*2 + end subroutine sub2 + + subroutine sub1(x,y) + integer, intent(in) :: x + integer, intent(out) :: y + + y = x + call sub2(y) + + + end subroutine sub1 + + +end module mod2 +! +! ! A simple generic linked list test program +! program list_test +! use linkedlist +! use data +! implicit none +! +! type(list_t), pointer :: ll => null() +! type(data_t), target :: dat_a +! type(data_t), target :: dat_b +! type(data_ptr) :: ptr +! +! ! Initialize two data objects +! dat_a%x = 17.5 +! dat_b%x = 3.0 +! +! ! Initialize the list with dat_a +! ptr%p => dat_a +! call list_init(ll, DATA=transfer(ptr, list_data)) +! print *, 'Initializing list with data:', ptr%p +! +! ! Insert dat_b into the list +! ptr%p => dat_b +! call list_insert(ll, DATA=transfer(ptr, list_data)) +! print *, 'Inserting node with data:', ptr%p +! +! ! Get the head node +! ptr = transfer(list_get(ll), ptr) +! print *, 'Head node data:', ptr%p +! +! ! Get the next node +! ptr = transfer(list_get(list_next(ll)), ptr) +! print *, 'Second node data:', ptr%p +! +! ! Free the list +! call list_free(ll) +! end program list_test + + + + + +program teste2 + use linkedlist + use data + use matprint + use mod3 + use mod1 + use mod2 + use mpi + + implicit none +! type(container), allocatable,dimension(:,:) :: m + integer,target :: d1,d2,d3,d4,d5,d6,d7,d8 + integer,pointer :: ptr + character(5) :: c + integer :: d + logical :: laux + real(dp) :: a(4), b(16) + real(dp), allocatable, dimension(:) :: r + integer :: i,j,jcell(10) + type(container), allocatable,dimension(:,:) :: m,n + integer ( kind = 4 ) id, root, dest + integer ( kind = 4 ) status(MPI_STATUS_SIZE) + integer ( kind = 4 ) ierr + integer ( kind = 4 ) p + integer ( kind = 4 ) tag + + allocate(m(2,1)) + allocate(n(2,1)) + d1 = 10 + ptr => d1 + call list_init(m(1,1)%list, DATA=transfer(ptr, list_data)) + + d2 =20 + ptr => d2 + call list_insert(m(1,1)%list, DATA=transfer(ptr, list_data)) + + d3 = 30 + ptr => d3 + call list_insert(m(1,1)%list, DATA=transfer(ptr, list_data)) + + d4 = 40 + ptr => d4 + call list_insert(m(1,1)%list, DATA=transfer(ptr, list_data)) + + d4 = 50 + ptr => d4 + call list_insert(m(1,1)%list, DATA=transfer(ptr, list_data)) + + + d5 = 1 + ptr => d5 + call list_init(m(2,1)%list, DATA=transfer(ptr, list_data)) + + d5 =2 + ptr => d5 + call list_insert(m(2,1)%list, DATA=transfer(ptr, list_data)) + + d5 = 3 + ptr => d5 + call list_insert(m(2,1)%list, DATA=transfer(ptr, list_data)) + + d5 = 4 + ptr => d5 + call list_insert(m(2,1)%list, DATA=transfer(ptr, list_data)) + + CALL RANDOM_SEED() + CALL RANDOM_NUMBER(r) + print*, r + + call MPI_Init ( ierr ) + call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr ) + call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr ) + + if (id == 0 .or. id == 1) then + a = [1,2,3,4] + ! CALL MPI_SEND(BUF, COUNT, DATATYPE, DEST, TAG, COMM, IERROR) + tag = 2 + dest = 2 + call MPI_send(a,4,MPI_DOUBLE_PRECISION,dest,tag,MPI_COMM_WORLD, ierr) + a = a*10 + tag = 3 + dest = 3 + call MPI_send(a,4,MPI_DOUBLE_PRECISION,dest,tag,MPI_COMM_WORLD, ierr) + + end if + + ! call MPI_Barrier(MPI_COMM_WORLD,ierr) + + ! if (id == 1 .or. id == 0) then + ! n(1,1) = m(1,1) + ! else + ! n(:,1) = m(:,1) + ! end if + ! deallocate(m) + + ! print*, ' ' + ! print*, '\nProcess',id + ! print*, ' ' + ! print*, 'ONE' + ! ptr = transfer(list_get(n(1,1)%list), ptr) + ! print *, 'Head node data:', ptr + + ! ! Get the next node + ! ptr = transfer(list_get(list_next(n(1,1)%list)), ptr) + ! print *, 'Second node data:', ptr + + ! print*, 'Associated n(1,2)?',associated(list_next(n(1,1)%list)) + ! ptr = transfer(list_get(list_next(n(1,1)%list)), ptr) + ! print *, ' data:', ptr + + call MPI_Finalize ( ierr ) + ! print*, 'OTHER' + ! ptr = transfer(list_get(m(2,1)%list), ptr) + ! print *, 'Head node data:', ptr + + ! ! Get the next node + ! ptr = transfer(list_get(list_next(m(2,1)%list)), ptr) + ! print *, 'Second node data:', ptr + + +! ptr => d5 +! print*, associated(m(2,1)%list) +! do d5 = 1,4 +! if (associated(m(2,1)%list)) then +! print*,d5 +! ! ptr => d5 +! print*, transfer(ptr, list_data) +! call list_insert(m(2,1)%list, DATA=transfer(ptr, list_data)) +! print*,d5 +! else +! print*,d5 +! ! ptr => d5 +! print*, transfer(ptr, list_data) +! call list_init(m(2,1)%list, DATA=transfer(ptr, list_data)) +! print*,d5 +! end if +! end do +! +! +! print*, "AAA" +! +! +! ! d = 3 +! ! ptr => d +! ! call list_init(m(2)%list, DATA=transfer(ptr, list_data)) +! ! +! ! d = 4 +! ! ptr => d +! ! call list_insert(m(2)%list, DATA=transfer(ptr, list_data)) +! +! +! +! ! ptr = transfer(list_get(m(2)%list), ptr) +! ! print *, 'Head node data:', ptr +! ! +! ! ! Get the next node +! ! ptr = transfer(list_get(list_next(m(2)%list)), ptr) +! ! print *, 'Second node data:', ptr +! ! +! print *, "OUTRA" +! +! ! ptr = transfer(list_get(m(2,1)%list), ptr) +! ! print *, 'Head node data:', ptr +! call testes(m,ptr,ptr) +! +! ptr = transfer(list_get(list_next(m(2,1)%list)), ptr) +! !print *, 'Second node data:', ptr +! +! node => list_next(m(2,1)%list) +! ! node2 => node +! ptr = transfer(list_get(node), ptr) +! print *, 'Second node data:', ptr +! +! node => list_next(node) +! node3 => list_next(node) +! +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Third node data:', ptr +! +! node => list_next(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fourth node data:', ptr +! +! node => list_next(node) +! if (associated(node)) then +! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fifth node data:', ptr +! end if +! +! print *, 'UMA' +! +! ptr = transfer(list_get(m(1,1)%list), ptr) +! print *, 'Head node data:', ptr +! +! ! ptr = transfer(list_get(list_next(m(1,1)%list)), ptr) +! ! print *, 'Second node data:', ptr +! +! node => list_next(m(1,1)%list) +! +! ptr = transfer(list_get(node), ptr) +! print *, 'Second node data:', ptr +! +! +! node2 => list_next(node) +! node => list_next(node) +! +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Third node data:', ptr +! +! node => list_next(node) +! +! ptr = transfer(list_get(node), ptr) +! print *, 'Fourth node data:', ptr +! +! node => list_next(node) +! print*,associated(node) +! if (associated(node)) then +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fifth node data:', ptr +! +! node => list_next(node) +! print*,associated(node) +! if (associated(node)) then +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Sixth node data:', ptr +! +! node => list_next(node) +! print*,associated(node) +! if (associated(node)) then +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Seventh node data:', ptr +! end if +! +! end if +! +! +! end if +! +! +! +! print *, "mudando o terceiro elemento" +! +! call testchange(node2,m) +! +! ptr = transfer(list_get(m(2,1)%list), ptr) +! print *, 'Head node data:', ptr +! +! ptr = transfer(list_get(list_next(m(2,1)%list)), ptr) +! !print *, 'Second node data:', ptr +! +! node => list_next(m(2,1)%list) +! ptr = transfer(list_get(node), ptr) +! print *, 'Second node data:', ptr +! +! node => list_next(node) +! +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Third node data:', ptr +! +! node => list_next(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fourth node data:', ptr +! +! node => list_next(node) +! if (associated(node)) then +! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fifth node data:', ptr +! end if +! +! print *, '\n original' +! +! ptr = transfer(list_get(m(1,1)%list), ptr) +! print *, 'Head node data:', ptr +! +! ptr = transfer(list_get(list_next(m(1,1)%list)), ptr) +! print *, 'Second node data:', ptr +! +! node => list_next(m(1,1)%list) +! ptr = transfer(list_get(node), ptr) +! ! print *, 'Second node data:', ptr +! +! node => list_next(node) +! +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Third node data:', ptr +! +! print*,'TESTE DOIDO' +! print*, ptr +! b = ptr +! print*, 'b=',b +! node => list_next(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fourth node data:', ptr +! +! node => list_next(node) +! print*,associated(node),'forth' +! if (associated(node)) then +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Fifth node data:', ptr +! +! node => list_next(node) +! print*,associated(node) +! if (associated(node)) then +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Sixth node data:', ptr +! +! node => list_next(node) +! print*,associated(node) +! if (associated(node)) then +! ! print*,associated(node) +! ptr = transfer(list_get(node), ptr) +! print *, 'Seventh node data:', ptr +! end if +! +! end if +! +! +! end if +! + + !ptr = transfer(list_get(node), ptr) + !print *, 'Fifth node data:', ptr + +! open(10,file='grupo1.csv',status="old") +! do i = 1,4 +! read(10,*) a(i,1),a(i,2) +! end do + + + +! call printm(a,'a') +! j = 2 +! jcell = (/((i*j),i = 1,10)/) +! +! laux = .true. +! do while (laux) +! print*, 'abc' +! end do +! PRINT *, "TESTE", (.NOT. 0) + c = 'abcde' + if ('abcde' == c) then + print*, "ola" + else + print*, "oi" + end if + + +end program teste2 +