c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c subroutine voigt_tabulation(eps_max,y_min,nint_y,file_out) implicit none include 'max.inc' include 'mpif.h' c c Purpose: to tabulate the Voigt profile as a function of x and y variables c integer strlen,idx integer nint_y integer i,j,ndiv integer n1(1:Nchunk_mx) integer n2(1:Nchunk_mx) integer it,Nit_max double precision f,pp double precision xres(1:Nppc_max,1:3) double precision xrest(1:Nppc_max,1:3,1:Nchunk_mx) double precision y_min,y_max,dx0 double precision eps_max character*(Nchar_mx) file_out double precision dp,dp_max,ys,y0 double precision y(1:Nchunk_mx),vP,lP,pmin,dy double precision pi integer a0,a25,a50,a75,a100 double precision apc c MPI integer ierr,np,pindex,nca integer stat(MPI_STATUS_SIZE) integer chunk,nchunks integer code0 integer proc,fin,fpf integer free(1:Nproc_mx) integer ch_idx(1:Nproc_mx) integer nsp c MPI 20 format(a3,e17.9,a9,i4) 50 format(100e17.9) CALL MPI_COMM_SIZE(MPI_COMM_WORLD,np,ierr) CALL MPI_COMM_RANK(MPI_COMM_WORLD,pindex,ierr) if (pindex.eq.0) then ! root process only pi=4.0D+0*datan(1.0D+0) ndiv=50 dx0=2.0D+0 Nit_max=100 pp=0.90D+0 c Look for y_max, the value of y at which the maximum relative difference between c the Lorentz profile and the Voigt profile is dp_max dp_max=1.0D-3 y0=30.0D+0 dy=2.0D+0 ys=y0 do it=1,Nit_max call Voigt(0.0D+0,y0,vP) lP=f(0.0D+0,y0)/dsqrt(pi) pmin=vP if (lP.lt.pmin) then pmin=lP endif dp=dabs(lP-vP)/pmin if (dp.lt.dp_max*pp) then ys=y0 y0=y0-dy endif if (dp.gt.dp_max) then y0=(y0+ys)/2.0D+0 dy=dy/2.0D+0 endif if ((dp.gt.dp_max*pp).and.(dp.lt.dp_max)) then y_max=y0 goto 126 endif enddo ! it 126 continue if (it.ge.Nit_max) then write(*,*) 'y_max could not be found' stop endif c Debug c write(*,*) 'y_max=',y_max c Debug nchunks=nint_y+1 do chunk=1,nchunks y(chunk)=y_min+(y_max-y_min)/(nchunks-1)*(chunk-1) enddo code0=5 ! data transfer do proc=1,np-1 call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) enddo c send integer data call MPI_BCAST(Nit_max,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ndiv,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) c send double precision data call MPI_BCAST(eps_max,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dx0,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(pp,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dp_max,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) do proc=1,np-1 free(proc)=1 enddo fin=0 chunk=0 a0=0 a25=0 a50=0 a75=0 a100=0 nca=0 do while(fin.eq.0) 444 continue call find_free_process(free,np-1,fpf,proc) if (fpf.eq.1) then ! a free child process found chunk=chunk+1 if (chunk.le.nchunks) then ch_idx(proc)=chunk code0=6 c Debug c write(*,*) 'Envoi de données à proc=',proc, c & chunk,'/',nchunks, c & ' [',numin(chunk),'-',numax(chunk),']' c Debug call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) call MPI_SEND(y(chunk),1,MPI_DOUBLE_PRECISION, & proc,1,MPI_COMM_WORLD,ierr) c Debug c write(*,*) 'chunk=',chunk,' y=',y(chunk),' running' c Debug free(proc)=0 ! child process index 'proc' is busy else ! chunk>nchunks free(proc)=2 ! means child process index 'proc' has been stopped endif goto 444 else ! no free child process found call stopped_processes(free,np-1,nsp) if (nsp.eq.np-1) then fin=1 goto 111 endif c wait for results call MPI_RECV(proc,1,MPI_INTEGER & ,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) idx=ch_idx(proc) nca=nca+1 c Debug c write(*,*) 'Réception de résultats depuis proc=',proc c Debug call MPI_RECV(n1(idx),1,MPI_INTEGER & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) call MPI_RECV(n2(idx),1,MPI_INTEGER & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) call MPI_RECV(xres,Nppc_max*3,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) apc=dble(nca)/dble(nchunks) call print_tabulation_progress(apc,a0,a25,a50,a75,a100) c Debug c write(*,*) 'nca:',nca,' /',nchunks,' apc=',apc c write(*,*) 'y=',y(idx),' n2=',n2(idx) c Debug free(proc)=1 ! child process index 'proc' is free again do i=1,n2(idx) do j=1,3 xrest(i,j,idx)=xres(i,j) enddo ! j enddo ! i endif ! fpf=0 or 1 111 continue enddo ! while fin=0 open(11,file=file_out(1:strlen(file_out))) do idx=1,nchunks write(11,20) 'y=',y(idx),' npoints=',n2(idx) do i=n1(idx)+1,n2(idx) write(11,50) (xrest(i,j,idx),j=1,3) enddo do i=1,n1(idx) write(11,50) (xrest(i,j,idx),j=1,3) enddo write(11,*) enddo ! idx close(11) write(*,*) file_out(1:strlen(file_out)), & ' was successfully generated' endif ! pindex return end subroutine print_tabulation_progress(apc,a0,a25,a50,a75,a100) implicit none double precision apc integer a0,a25,a50,a75,a100 if ((apc.ge.0.0D+0).and.(apc.lt.0.25D+0).and.(a0.eq.0)) & then write(*,*) 'achieved: 0%' a0=1 endif if ((apc.ge.0.25D+0).and.(apc.lt.0.50D+0).and.(a25.eq.0)) & then write(*,*) 'achieved: 25%' a25=1 endif if ((apc.ge.0.50D+0).and.(apc.lt.0.75D+0).and.(a50.eq.0)) & then write(*,*) 'achieved: 50%' a50=1 endif if ((apc.ge.0.75D+0).and.(apc.lt.1.00D+0).and.(a75.eq.0)) & then write(*,*) 'achieved: 75%' a75=1 endif if ((apc.ge.1.0D+0).and.(a100.eq.0)) & then write(*,*) 'achieved: 100%' a100=1 endif return end subroutine voigt_tabulation_1y(y,eps_max,dx0,Nit_max,pp,ndiv, & dp_max,n1,n2,xres) implicit none include 'max.inc' integer Nit_max,ndiv integer n1,n2,i,fin,it double precision y,eps_max,dx0,pp double precision x1p,x1,emax,epsm,x,xs,dx,x2,f double precision xres(1:Nppc_max,1:3) double precision vP,lP,pi,pmin,dp,dp_max pi=4.0D+0*datan(1.0D+0) n1=0 i=0 x1p=0.0D+0 fin=0 do while (fin.eq.0) x1=x1p c set variable emax if (dabs(x1).le.1.0D+0) then emax=eps_max/1.5D+1 else emax=eps_max endif epsm=0.0D+0 it=0 x=x1+dx0 xs=x1 dx=dx0 do it=1,Nit_max call find_epsm_vt(x1,x,y,epsm) if (epsm.lt.emax*pp) then xs=x x=x+dx dx=dx*1.1D+0 endif if (epsm.gt.emax) then x=(xs+x)/2.0D+0 dx=dx/1.1D+0 endif if ((epsm.gt.emax*pp).and.(epsm.lt.emax)) then x2=x goto 123 endif enddo ! it if (it.ge.Nit_max) then write(*,*) 'error: Nit_max reached' stop endif 123 continue x1p=x1+(x2-x1)/ndiv i=i+1 if (i.ge.Nppc_max) then write(*,*) 'Error: Nppc_max reached' stop endif xres(i,1)=x1 xres(i,2)=x1p xres(i,3)=x2 dx0=1.1D+0*(x2-x1) call Voigt(x1p,y,vP) lP=f(x1p,y)/dsqrt(pi) pmin=vP if (lP.lt.pmin) then pmin=lP endif dp=dabs(lP-vP)/pmin if (dp.lt.dp_max) then fin=1 endif enddo n1=i x1p=-xres(n1,2) fin=0 do while (fin.eq.0) x1=x1p c set variable emax if (dabs(x1).le.1.0D+0) then emax=eps_max/1.5D+1 else emax=eps_max endif epsm=0.0D+0 it=0 x=x1+dx0 xs=x1 dx=dx0 do it=1,Nit_max call find_epsm_vt(x1,x,y,epsm) if (epsm.lt.emax*pp) then xs=x x=x+dx endif if (epsm.gt.emax) then x=(xs+x)/2.0D+0 dx=dx/2.0D+0 endif if ((epsm.gt.emax*pp).and.(epsm.lt.emax)) then x2=x goto 124 endif enddo ! it if (it.ge.Nit_max) then write(*,*) 'error: Nit_max reached' stop endif 124 continue x1p=x1+(x2-x1)/ndiv dx0=1.1D+0*(x2-x1) if ((x1.lt.0.0D+0) & .and.(x2.gt.0.0D+0) & .and.(epsm.lt.emax)) then x2=0.0D+0 if (x1p.gt.0.0D+0) then x1p=0.0D+0 fin=1 endif endif i=i+1 if (i.ge.Nppc_max) then write(*,*) 'Error: Nppc_max reached' stop endif xres(i,1)=x1 xres(i,2)=x1p xres(i,3)=x2 enddo n2=i return end subroutine find_epsm_vt(x1,x,y,epsm) implicit none include 'max.inc' double precision fs,vP,min double precision x,y,x1,epsm,xl double precision eps(1:Nmax_vt) integer N,i,imax c Debug if (x.lt.x1) then write(*,*) 'This is routine find_epsm_vt' write(*,*) 'x1=',x1,' x=',x endif c Debug N=nint((x-x1)*1.0D+4) if (N.gt.Nmax_vt) then do while (N.gt.Nmax_vt) N=N/10 enddo endif do i=1,N xl=x1+(x-x1)/(N-1)*(i-1) call fstar(x1,x,xl,y,fs) call Voigt(xl,y,vP) min=fs if (vP.lt.min) then min=vP endif eps(i)=dabs(vP-fs)/min enddo call find_max_vt(eps,N,epsm,imax) return end double precision function f(x,y) implicit none double precision x,y f=y/(y**2+x**2) return end subroutine fstar(x1,x2,x,y,fs) implicit none double precision x1,x2,x,y,fs,f1,f2 call Voigt(x1,y,f1) call Voigt(x2,y,f2) fs=f1+(x-x1)*(f2-f1)/(x2-x1) return end