From a38db3320feb55f87333840f852c90df7c0f4573 Mon Sep 17 00:00:00 2001 From: Pierre Navaro Date: Mon, 25 Mar 2024 13:25:11 +0100 Subject: [PATCH] run fprettify and JuliaFormatter --- commun.F90 | 79 +++++---- fdtd.jl | 48 +++--- maxwell_mpi.F90 | 364 +++++++++++++++++++-------------------- maxwell_mpi.jl | 233 ++++++++++++++----------- maxwell_serial.F90 | 133 +++++++------- maxwell_serial.jl | 98 ++++++----- mesh.jl | 55 +++--- plot_mpi.jl | 61 ++++--- read_data.jl | 49 ------ ring.jl | 38 ++-- ring_with_distributed.jl | 40 ++--- solveur_yee.F90 | 100 ++++++----- sorties.F90 | 124 ++++++------- 13 files changed, 695 insertions(+), 727 deletions(-) delete mode 100644 read_data.jl diff --git a/commun.F90 b/commun.F90 index 803d138..6274f64 100644 --- a/commun.F90 +++ b/commun.F90 @@ -1,60 +1,59 @@ module commun -type mesh_fields - real(8), dimension(:,:), pointer :: ex, ey - real(8), dimension(:,:), pointer :: bz -end type mesh_fields + type mesh_fields + real(8), dimension(:, :), pointer :: ex, ey + real(8), dimension(:, :), pointer :: bz + end type mesh_fields -real(8) :: c, csq -real(8) :: pi + real(8) :: c, csq + real(8) :: pi -integer :: md, nd -integer :: nx, ny -integer :: nstep, nstepmax -integer :: idiag + integer :: md, nd + integer :: nx, ny + integer :: nstep, nstepmax + integer :: idiag -integer, private :: i, j + integer, private :: i, j -real(8) :: dt, dx, dy -real(8) :: dimx, dimy -real(8) :: cfl -real(8) :: tfinal -real(8) :: omega + real(8) :: dt, dx, dy + real(8) :: dimx, dimy + real(8) :: cfl + real(8) :: tfinal + real(8) :: omega contains -subroutine readin( ) + subroutine readin() -implicit none + implicit none -namelist/donnees/ cfl, & !nbre de Courant - tfinal, & !duree maxi - nstepmax, & !nbre d'iterations maxi - idiag, & !frequence des diagnostics - md,nd, & !nombre d'onde de la solution initiale - nx, ny, & !nombre de points dans les deux directions - dimx, dimy !dimensions du domaine + namelist /donnees/ cfl, & !nbre de Courant + tfinal, & !duree maxi + nstepmax, & !nbre d'iterations maxi + idiag, & !frequence des diagnostics + md, nd, & !nombre d'onde de la solution initiale + nx, ny, & !nombre de points dans les deux directions + dimx, dimy !dimensions du domaine !***Initialisation des valeurs pas default -pi = 4.*atan(1.) + pi = 4.*atan(1.) -c = 1d0 !celerite de la lumiere -csq = c * c + c = 1d0 !celerite de la lumiere + csq = c*c -open(10,file="input_data",status='old') -read(10,donnees) -close(10) + open (10, file="input_data", status='old') + read (10, donnees) + close (10) -write(*,"(a,g12.3)")" largeur dimx = ", dimx -write(*,"(a,g12.3)")" longueur dimy = ", dimy -write(*,"(a,g12.3)")" temps = ", tfinal -write(*,"(a,g12.3)")" nombre nx = ", nx -write(*,"(a,g12.3)")" nombre ny = ", ny -write(*,"(a,g12.3)")" nombre de Courant = ", cfl -write(*,"(a,g12.3)")" frequence des diagnostics = ", idiag + write (*, "(a,g12.3)") " largeur dimx = ", dimx + write (*, "(a,g12.3)") " longueur dimy = ", dimy + write (*, "(a,g12.3)") " temps = ", tfinal + write (*, "(a,g12.3)") " nombre nx = ", nx + write (*, "(a,g12.3)") " nombre ny = ", ny + write (*, "(a,g12.3)") " nombre de Courant = ", cfl + write (*, "(a,g12.3)") " frequence des diagnostics = ", idiag - -end subroutine readin + end subroutine readin end module commun diff --git a/fdtd.jl b/fdtd.jl index db77d44..48f7f79 100644 --- a/fdtd.jl +++ b/fdtd.jl @@ -1,46 +1,46 @@ -function faraday!( fields, ix, jx, iy, jy, dt ) +function faraday!(fields, ix, jx, iy, jy, dt) dx, dy = fields.mesh.dx, fields.mesh.dy - for j=iy:jy, i=ix:jx - dex_dy = (fields.ex[i,j+1]-fields.ex[i,j]) / dy - dey_dx = (fields.ey[i+1,j]-fields.ey[i,j]) / dx - fields.bz[i,j] = fields.bz[i,j] + dt * (dex_dy - dey_dx) + for j = iy:jy, i = ix:jx + dex_dy = (fields.ex[i, j+1] - fields.ex[i, j]) / dy + dey_dx = (fields.ey[i+1, j] - fields.ey[i, j]) / dx + fields.bz[i, j] = fields.bz[i, j] + dt * (dex_dy - dey_dx) end end -function ampere_maxwell!( fields, ix, jx, iy, jy, dt ) +function ampere_maxwell!(fields, ix, jx, iy, jy, dt) dx, dy = fields.mesh.dx, fields.mesh.dy - for j=iy+1:jy, i=ix:jx - dbz_dy = (fields.bz[i,j]-fields.bz[i,j-1]) / dy - fields.ex[i,j] = fields.ex[i,j] + dt*csq*dbz_dy + for j = iy+1:jy, i = ix:jx + dbz_dy = (fields.bz[i, j] - fields.bz[i, j-1]) / dy + fields.ex[i, j] = fields.ex[i, j] + dt * csq * dbz_dy end - for j=iy:jy, i=ix+1:jx - dbz_dx = (fields.bz[i,j]-fields.bz[i-1,j]) / dx - fields.ey[i,j] = fields.ey[i,j] - dt*csq*dbz_dx + for j = iy:jy, i = ix+1:jx + dbz_dx = (fields.bz[i, j] - fields.bz[i-1, j]) / dx + fields.ey[i, j] = fields.ey[i, j] - dt * csq * dbz_dx end -end +end function periodic_bc!(fields, ix, jx, iy, jy, dt) for i = ix:jx - dbz_dy = (fields.bz[i,iy]-fields.bz[i,jy]) / fields.mesh.dy - fields.ex[i,iy] = fields.ex[i,iy] + dt*csq*dbz_dy - fields.ex[i,jy+1] = fields.ex[i,iy] - end - + dbz_dy = (fields.bz[i, iy] - fields.bz[i, jy]) / fields.mesh.dy + fields.ex[i, iy] = fields.ex[i, iy] + dt * csq * dbz_dy + fields.ex[i, jy+1] = fields.ex[i, iy] + end + for j = iy:jy - dbz_dx = (fields.bz[ix,j]-fields.bz[jx,j]) / fields.mesh.dx - fields.ey[ix,j] = fields.ey[ix,j] - dt*csq*dbz_dx - fields.ey[jx+1,j] = fields.ey[ix,j] + dbz_dx = (fields.bz[ix, j] - fields.bz[jx, j]) / fields.mesh.dx + fields.ey[ix, j] = fields.ey[ix, j] - dt * csq * dbz_dx + fields.ey[jx+1, j] = fields.ey[ix, j] end - fields.bz[:,jy+1] .= fields.bz[:,iy] - fields.bz[jx+1,:] .= fields.bz[ix,:] - + fields.bz[:, jy+1] .= fields.bz[:, iy] + fields.bz[jx+1, :] .= fields.bz[ix, :] + end diff --git a/maxwell_mpi.F90 b/maxwell_mpi.F90 index bbe5e64..7506578 100644 --- a/maxwell_mpi.F90 +++ b/maxwell_mpi.F90 @@ -1,226 +1,222 @@ program maxyee_mpi -use mpi -use commun -use sorties - -implicit none - -type(mesh_fields) :: fields -integer :: i, j, iproc, istep, iplot -real(8) :: tcpu1, tcpu2, time, err_l2, sum_l2, xp, yp -real(8) :: x, y, dtloc -real(8) :: th_bz -integer,dimension(MPI_STATUS_SIZE) :: statut -integer :: rang, nproc, code, comm2d -integer,parameter :: tag=1111 -integer,dimension(8) :: voisin -integer,parameter :: N =1, S =2, W =3, E =4 -integer,parameter :: ndims = 2 -integer,dimension(ndims) :: dims, coords -logical :: reorder -logical,dimension(ndims) :: periods -integer :: nxp, nyp, mx, my -real(8) :: dex_dx, dey_dy -real(8) :: dex_dy, dey_dx -real(8) :: dbz_dx, dbz_dy -integer :: type_ligne, type_colonne + use mpi + use commun + use sorties + + implicit none + + type(mesh_fields) :: fields + integer :: i, j, iproc, istep, iplot + real(8) :: tcpu1, tcpu2, time, err_l2, sum_l2, xp, yp + real(8) :: x, y, dtloc + real(8) :: th_bz + integer, dimension(MPI_STATUS_SIZE) :: statut + integer :: rang, nproc, code, comm2d + integer, parameter :: tag = 1111 + integer, dimension(8) :: voisin + integer, parameter :: N = 1, S = 2, W = 3, E = 4 + integer, parameter :: ndims = 2 + integer, dimension(ndims) :: dims, coords + logical :: reorder + logical, dimension(ndims) :: periods + integer :: nxp, nyp, mx, my + real(8) :: dex_dx, dey_dy + real(8) :: dex_dy, dey_dx + real(8) :: dbz_dx, dbz_dy + integer :: type_ligne, type_colonne !Initialisation de MPI -CALL MPI_INIT(code) -CALL MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,code) -CALL MPI_COMM_RANK(MPI_COMM_WORLD,rang,code) -CALL MPI_BARRIER(MPI_COMM_WORLD, code) + CALL MPI_INIT(code) + CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, code) + CALL MPI_COMM_RANK(MPI_COMM_WORLD, rang, code) + CALL MPI_BARRIER(MPI_COMM_WORLD, code) -tcpu1 = MPI_WTIME() + tcpu1 = MPI_WTIME() !Nombre de processus suivant x et y -dims(:) = 0 -CALL MPI_DIMS_CREATE(nproc,ndims,dims,code) -nxp = dims(1) -nyp = dims(2) + dims(:) = 0 + CALL MPI_DIMS_CREATE(nproc, ndims, dims, code) + nxp = dims(1) + nyp = dims(2) !Creation de la grille 2D periodique en x et y -periods(1) = .true. -periods(2) = .true. -reorder = .true. + periods(1) = .true. + periods(2) = .true. + reorder = .true. -CALL MPI_CART_CREATE(MPI_COMM_WORLD,ndims,dims,periods,reorder,comm2d,code) + CALL MPI_CART_CREATE(MPI_COMM_WORLD, ndims, dims, periods, reorder, comm2d, code) -!Initialisation du tableau voisin -voisin(:) = MPI_PROC_NULL +!Initialisation du tableau voisin + voisin(:) = MPI_PROC_NULL !Recherche de mes voisins Sud et Nord -CALL MPI_CART_SHIFT(comm2d,0,1,voisin(N),voisin(S),code) + CALL MPI_CART_SHIFT(comm2d, 0, 1, voisin(N), voisin(S), code) !Recherche de mes voisins Ouest et Est -CALL MPI_CART_SHIFT(comm2d,1,1,voisin(W),voisin(E),code) + CALL MPI_CART_SHIFT(comm2d, 1, 1, voisin(W), voisin(E), code) !Connaitre mes coordonnees dans la topologie -CALL MPI_COMM_RANK(comm2d,rang,code) -CALL MPI_CART_COORDS(comm2d,rang,ndims,coords,code) - -if (rang == 0) call readin( ) - -call MPI_BCAST(nx, 1,MPI_INTEGER, 0,comm2d,code) -call MPI_BCAST(ny, 1,MPI_INTEGER, 0,comm2d,code) -call MPI_BCAST(nd, 1,MPI_INTEGER, 0,comm2d,code) -call MPI_BCAST(md, 1,MPI_INTEGER, 0,comm2d,code) -call MPI_BCAST(idiag, 1,MPI_INTEGER, 0,comm2d,code) -call MPI_BCAST(dimx, 1,MPI_REAL8, 0,comm2d,code) -call MPI_BCAST(dimy, 1,MPI_REAL8, 0,comm2d,code) -call MPI_BCAST(pi, 1,MPI_REAL8, 0,comm2d,code) -call MPI_BCAST(csq, 1,MPI_REAL8, 0,comm2d,code) -call MPI_BCAST(c, 1,MPI_REAL8, 0,comm2d,code) - -mx = nx/nxp; dx = dimx/nx -my = ny/nyp; dy = dimy/ny - -call MPI_BCAST(cfl, 1,MPI_REAL8, 0,comm2d,code) -dtloc = cfl / sqrt (1./(dx*dx)+1./(dy*dy)) / c -call MPI_ALLREDUCE(dtloc,dt,1,MPI_REAL8,MPI_MAX,comm2d,code) - -call MPI_BCAST(tfinal,1,MPI_REAL8, 0,comm2d,code) -call MPI_BCAST(nstepmax,1,MPI_INTEGER, 0,comm2d,code) -nstep = floor(tfinal/dt) -if( nstep > nstepmax ) nstep = nstepmax - -allocate(fields%ex(mx+1,my+1)); fields%ex(:,:) = 0.0d0 -allocate(fields%ey(mx+1,my+1)); fields%ey(:,:) = 0.0d0 -allocate(fields%bz(mx+1,my+1)); fields%bz(:,:) = 0.0d0 - -xp = coords(1) * dimx / nxp -yp = coords(2) * dimy / nyp - -print*, "proc= ",rang,": ",mx,"x",my," mes coords sont ",coords(:) -print*, "proc= ",rang,": (xp,yp) = (",sngl(xp),";", sngl(yp), ")" -print*, "proc= ",rang,": mes voisins sont ", voisin(1:4) -print*, "proc= ",rang,": Nombre d'iteration nstep = ", nstep -print*, "proc= ",rang,": dx = ", sngl(dx), " dy = ", sngl(dy), " dt = ", sngl(dt) - + CALL MPI_COMM_RANK(comm2d, rang, code) + CALL MPI_CART_COORDS(comm2d, rang, ndims, coords, code) + + if (rang == 0) call readin() + + call MPI_BCAST(nx, 1, MPI_INTEGER, 0, comm2d, code) + call MPI_BCAST(ny, 1, MPI_INTEGER, 0, comm2d, code) + call MPI_BCAST(nd, 1, MPI_INTEGER, 0, comm2d, code) + call MPI_BCAST(md, 1, MPI_INTEGER, 0, comm2d, code) + call MPI_BCAST(idiag, 1, MPI_INTEGER, 0, comm2d, code) + call MPI_BCAST(dimx, 1, MPI_REAL8, 0, comm2d, code) + call MPI_BCAST(dimy, 1, MPI_REAL8, 0, comm2d, code) + call MPI_BCAST(pi, 1, MPI_REAL8, 0, comm2d, code) + call MPI_BCAST(csq, 1, MPI_REAL8, 0, comm2d, code) + call MPI_BCAST(c, 1, MPI_REAL8, 0, comm2d, code) + + mx = nx/nxp; dx = dimx/nx + my = ny/nyp; dy = dimy/ny + + call MPI_BCAST(cfl, 1, MPI_REAL8, 0, comm2d, code) + dtloc = cfl/sqrt(1./(dx*dx) + 1./(dy*dy))/c + call MPI_ALLREDUCE(dtloc, dt, 1, MPI_REAL8, MPI_MAX, comm2d, code) + + call MPI_BCAST(tfinal, 1, MPI_REAL8, 0, comm2d, code) + call MPI_BCAST(nstepmax, 1, MPI_INTEGER, 0, comm2d, code) + nstep = floor(tfinal/dt) + if (nstep > nstepmax) nstep = nstepmax + + allocate (fields%ex(mx + 1, my + 1)); fields%ex(:, :) = 0.0d0 + allocate (fields%ey(mx + 1, my + 1)); fields%ey(:, :) = 0.0d0 + allocate (fields%bz(mx + 1, my + 1)); fields%bz(:, :) = 0.0d0 + + xp = coords(1)*dimx/nxp + yp = coords(2)*dimy/nyp + + print *, "proc= ", rang, ": ", mx, "x", my, " mes coords sont ", coords(:) + print *, "proc= ", rang, ": (xp,yp) = (", sngl(xp), ";", sngl(yp), ")" + print *, "proc= ", rang, ": mes voisins sont ", voisin(1:4) + print *, "proc= ", rang, ": Nombre d'iteration nstep = ", nstep + print *, "proc= ", rang, ": dx = ", sngl(dx), " dy = ", sngl(dy), " dt = ", sngl(dt) !type colonne -CALL MPI_TYPE_CONTIGUOUS(mx+1,MPI_REAL8,type_colonne,code) -CALL MPI_TYPE_COMMIT(type_colonne,code) + CALL MPI_TYPE_CONTIGUOUS(mx + 1, MPI_REAL8, type_colonne, code) + CALL MPI_TYPE_COMMIT(type_colonne, code) !type ligne -CALL MPI_TYPE_VECTOR(my+1,1,mx+1,MPI_REAL8,type_ligne,code) -CALL MPI_TYPE_COMMIT(type_ligne,code) + CALL MPI_TYPE_VECTOR(my + 1, 1, mx + 1, MPI_REAL8, type_ligne, code) + CALL MPI_TYPE_COMMIT(type_ligne, code) + CALL FLUSH (6) + call MPI_BARRIER(MPI_COMM_WORLD, code) -CALL FLUSH(6) -call MPI_BARRIER(MPI_COMM_WORLD, code) + time = 0. + iplot = 0 -time = 0. -iplot = 0 +!Initialisation des champs + fields%ex(:, :) = 0d0; fields%ey(:, :) = 0d0; fields%bz(:, :) = 0d0 -!Initialisation des champs -fields%ex(:,:) = 0d0; fields%ey(:,:) = 0d0; fields%bz(:,:) = 0d0 + xp = coords(1)*dimx/nxp + yp = coords(2)*dimy/nyp -xp = coords(1) * dimx /nxp -yp = coords(2) * dimy /nyp - -omega = c * sqrt((md*pi/dimx)**2+(nd*pi/dimy)**2) -do j=1,my - do i=1,mx - x = xp + (i-.5) * dx - y = yp + (j-.5) * dy - fields%bz(i,j) = - cos(md*pi*x/dimx) & - * cos(nd*pi*y/dimy) & - * cos(omega*(-0.5*dt)) - end do -end do + omega = c*sqrt((md*pi/dimx)**2 + (nd*pi/dimy)**2) + do j = 1, my + do i = 1, mx + x = xp + (i - .5)*dx + y = yp + (j - .5)*dy + fields%bz(i, j) = -cos(md*pi*x/dimx) & + *cos(nd*pi*y/dimy) & + *cos(omega*(-0.5*dt)) + end do + end do -do istep = 1, nstep !*** Loop over time + do istep = 1, nstep !*** Loop over time - !E(n) [1:mx]*[1:my] --> B(n+1/2) [1:mx-1]*[1:my-1] + !E(n) [1:mx]*[1:my] --> B(n+1/2) [1:mx-1]*[1:my-1] - do j=1,my - do i=1,mx - dex_dy = (fields%ex(i,j+1)-fields%ex(i,j)) / dy - dey_dx = (fields%ey(i+1,j)-fields%ey(i,j)) / dx - fields%bz(i,j) = fields%bz(i,j) + dt * (dex_dy - dey_dx) - end do - end do + do j = 1, my + do i = 1, mx + dex_dy = (fields%ex(i, j + 1) - fields%ex(i, j))/dy + dey_dx = (fields%ey(i + 1, j) - fields%ey(i, j))/dx + fields%bz(i, j) = fields%bz(i, j) + dt*(dex_dy - dey_dx) + end do + end do - time = time + 0.5*dt + time = time + 0.5*dt - CALL MPI_SENDRECV(fields%bz( 1, 1),1,type_ligne,voisin(N),tag, & - fields%bz(mx+1, 1),1,type_ligne,voisin(S),tag, & - comm2d, statut, code) + CALL MPI_SENDRECV(fields%bz(1, 1), 1, type_ligne, voisin(N), tag, & + fields%bz(mx + 1, 1), 1, type_ligne, voisin(S), tag, & + comm2d, statut, code) - CALL MPI_SENDRECV(fields%bz( 1, 1),1,type_colonne,voisin(W),tag, & - fields%bz( 1,my+1),1,type_colonne,voisin(E),tag, & - comm2d, statut, code) + CALL MPI_SENDRECV(fields%bz(1, 1), 1, type_colonne, voisin(W), tag, & + fields%bz(1, my + 1), 1, type_colonne, voisin(E), tag, & + comm2d, statut, code) + !Bz(n+1/2) [1:mx]*[1:my] --> Ex(n+1) [1:mx]*[2:my] + !Bz(n+1/2) [1:mx]*[1:my] --> Ey(n+1) [2:mx]*[1:my] + do j = 2, my + 1 + do i = 1, mx + dbz_dy = (fields%bz(i, j) - fields%bz(i, j - 1))/dy + fields%ex(i, j) = fields%ex(i, j) + dt*csq*dbz_dy + end do + end do - !Bz(n+1/2) [1:mx]*[1:my] --> Ex(n+1) [1:mx]*[2:my] - !Bz(n+1/2) [1:mx]*[1:my] --> Ey(n+1) [2:mx]*[1:my] - do j=2,my+1 - do i=1,mx - dbz_dy = (fields%bz(i,j)-fields%bz(i,j-1)) / dy - fields%ex(i,j) = fields%ex(i,j) + dt*csq*dbz_dy - end do - end do - - do j=1,my - do i=2,mx+1 - dbz_dx = (fields%bz(i,j)-fields%bz(i-1,j)) / dx - fields%ey(i,j) = fields%ey(i,j) - dt*csq*dbz_dx - end do - end do + do j = 1, my + do i = 2, mx + 1 + dbz_dx = (fields%bz(i, j) - fields%bz(i - 1, j))/dx + fields%ey(i, j) = fields%ey(i, j) - dt*csq*dbz_dx + end do + end do - !Envoi au voisin E et reception du voisin W + !Envoi au voisin E et reception du voisin W - CALL MPI_SENDRECV(fields%ex( 1,my+1),1,type_colonne,voisin(E),tag, & - fields%ex( 1, 1),1,type_colonne,voisin(W),tag, & - comm2d, statut, code) + CALL MPI_SENDRECV(fields%ex(1, my + 1), 1, type_colonne, voisin(E), tag, & + fields%ex(1, 1), 1, type_colonne, voisin(W), tag, & + comm2d, statut, code) !!Envoi au voisin S et reception du voisin N - CALL MPI_SENDRECV(fields%ey(mx+1, 1),1,type_ligne,voisin(S),tag, & - fields%ey( 1, 1),1,type_ligne,voisin(N),tag, & - comm2d, statut, code) - - !-----------------------------------------------------------! - ! Sorties graphiques (les champs sont connus au temps n) ! - ! pour le solveur de MAXWELL ! - !-----------------------------------------------------------! - err_l2 = 0.0 - do j = 1, my - do i = 1, mx - x = xp + (i-.5) * dx - y = yp + (j-.5) * dy - th_bz = - cos(md*pi*x/dimx) & - * cos(nd*pi*y/dimy) & - * cos(omega*time) - err_l2 = err_l2 + (fields%bz(i,j) - th_bz)**2 - end do - end do - - - if ( mod(istep,idiag) == 0.0) then - iplot = iplot + 1 - call plot_fields(rang, nproc, fields, 1, mx, 1, my, xp, yp, iplot, time ) + CALL MPI_SENDRECV(fields%ey(mx + 1, 1), 1, type_ligne, voisin(S), tag, & + fields%ey(1, 1), 1, type_ligne, voisin(N), tag, & + comm2d, statut, code) + + !-----------------------------------------------------------! + ! Sorties graphiques (les champs sont connus au temps n) ! + ! pour le solveur de MAXWELL ! + !-----------------------------------------------------------! + err_l2 = 0.0 + do j = 1, my + do i = 1, mx + x = xp + (i - .5)*dx + y = yp + (j - .5)*dy + th_bz = -cos(md*pi*x/dimx) & + *cos(nd*pi*y/dimy) & + *cos(omega*time) + err_l2 = err_l2 + (fields%bz(i, j) - th_bz)**2 + end do + end do + + if (mod(istep, idiag) == 0.0) then + iplot = iplot + 1 + call plot_fields(rang, nproc, fields, 1, mx, 1, my, xp, yp, iplot, time) + end if + + time = time + 0.5*dt + + end do ! next time step + + call MPI_REDUCE(err_l2, sum_l2, 1, MPI_REAL8, MPI_SUM, 0, comm2d, code) + + if (rang == 0) then + write (*, "(10x,' istep = ',I6)", advance="no") istep + write (*, "(' time = ',e17.3,' ns')", advance="no") time + write (*, "(' erreur L2 = ',g10.3)") sqrt(sum_l2) end if - time = time + 0.5*dt - -end do ! next time step - -call MPI_REDUCE (err_l2,sum_l2,1,MPI_REAL8,MPI_SUM,0,comm2d,code) - -if (rang == 0) then - write(*,"(10x,' istep = ',I6)",advance="no") istep - write(*,"(' time = ',e17.3,' ns')",advance="no") time - write(*,"(' erreur L2 = ',g10.3)") sqrt(sum_l2) -end if - -call FLUSH(6) -call MPI_BARRIER(MPI_COMM_WORLD, code) -tcpu2 = MPI_WTIME() -if (rang == 0) & - write(*,"(//10x,' Temps CPU = ', G15.3, ' sec' )") (tcpu2-tcpu1)*nproc + call FLUSH (6) + call MPI_BARRIER(MPI_COMM_WORLD, code) + tcpu2 = MPI_WTIME() + if (rang == 0) & + write (*, "(//10x,' Temps CPU = ', G15.3, ' sec' )") (tcpu2 - tcpu1)*nproc !Desactivation de MPI -CALL MPI_FINALIZE(code) -stop + CALL MPI_FINALIZE(code) + stop 1000 format(10f8.4) diff --git a/maxwell_mpi.jl b/maxwell_mpi.jl index 703c6ea..60e924a 100644 --- a/maxwell_mpi.jl +++ b/maxwell_mpi.jl @@ -6,69 +6,69 @@ const csq = c * c struct Mesh - nx :: Int64 - dx :: Float64 - ny :: Int64 - dy :: Float64 + nx::Int64 + dx::Float64 + ny::Int64 + dy::Float64 end struct MeshFields - mesh :: Mesh - ex :: Array{Float64, 2} - ey :: Array{Float64, 2} - bz :: Array{Float64, 2} + mesh::Mesh + ex::Array{Float64,2} + ey::Array{Float64,2} + bz::Array{Float64,2} - function MeshFields( mesh ) + function MeshFields(mesh) nx, ny = mesh.nx, mesh.ny - ex = zeros(Float64, (nx,ny+1)) - ey = zeros(Float64, (nx+1,ny)) - bz = zeros(Float64, (nx+1,ny+1)) - new( mesh, ex, ey, bz) + ex = zeros(Float64, (nx, ny + 1)) + ey = zeros(Float64, (nx + 1, ny)) + bz = zeros(Float64, (nx + 1, ny + 1)) + new(mesh, ex, ey, bz) end -end +end -function faraday!( fields, dt ) +function faraday!(fields, dt) dx, dy = fields.mesh.dx, fields.mesh.dy nx, ny = fields.mesh.nx, fields.mesh.ny - for j=1:ny, i=1:nx - dex_dy = (fields.ex[i,j+1]-fields.ex[i,j]) / dy - dey_dx = (fields.ey[i+1,j]-fields.ey[i,j]) / dx - fields.bz[i,j] = fields.bz[i,j] + dt * (dex_dy - dey_dx) + for j = 1:ny, i = 1:nx + dex_dy = (fields.ex[i, j+1] - fields.ex[i, j]) / dy + dey_dx = (fields.ey[i+1, j] - fields.ey[i, j]) / dx + fields.bz[i, j] = fields.bz[i, j] + dt * (dex_dy - dey_dx) end end -function ampere_maxwell!( fields, dt ) +function ampere_maxwell!(fields, dt) dx, dy = fields.mesh.dx, fields.mesh.dy nx, ny = fields.mesh.nx, fields.mesh.ny - for j=2:ny+1, i=1:nx - dbz_dy = (fields.bz[i,j]-fields.bz[i,j-1]) / dy - fields.ex[i,j] = fields.ex[i,j] + dt*csq*dbz_dy + for j = 2:ny+1, i = 1:nx + dbz_dy = (fields.bz[i, j] - fields.bz[i, j-1]) / dy + fields.ex[i, j] = fields.ex[i, j] + dt * csq * dbz_dy end - for j=1:ny, i=2:nx+1 - dbz_dx = (fields.bz[i,j]-fields.bz[i-1,j]) / dx - fields.ey[i,j] = fields.ey[i,j] - dt*csq*dbz_dx + for j = 1:ny, i = 2:nx+1 + dbz_dx = (fields.bz[i, j] - fields.bz[i-1, j]) / dx + fields.ey[i, j] = fields.ey[i, j] - dt * csq * dbz_dx end -end +end include("plot_mpi.jl") -function main( nstep ) +function main(nstep) - cfl = 0.4 # Courant-Friedrich-Levy - tfinal = 1. # final time + cfl = 0.4 # Courant-Friedrich-Levy + tfinal = 1.0 # final time nstepmax = 1000 # max steps md = 2 # md : wave number x (initial condition) nd = 2 # nd : wave number y (initial condition) @@ -84,108 +84,133 @@ function main( nstep ) rank = MPI.Comm_rank(comm) MPI.Barrier(comm) - + tcpu1 = MPI.Wtime() - + dims = [0 0] ndims = length(dims) periods = [1 1] reorder = 1 - + MPI.Dims_create!(proc, dims) comm2d = MPI.Cart_create(comm, dims, periods, reorder) - + @assert MPI.Comm_size(comm2d) == proc - - north, south = MPI.Cart_shift(comm2d,0,1) - west, east = MPI.Cart_shift(comm2d,1,1) - + + north, south = MPI.Cart_shift(comm2d, 0, 1) + west, east = MPI.Cart_shift(comm2d, 1, 1) + coords = MPI.Cart_coords(comm2d) - + nxp, nyp = dims mx = nx ÷ nxp my = ny ÷ nyp - dt = cfl / sqrt(1/dx^2+1/dy^2) / c - - nstep = min( nstepmax, nstep) - + dt = cfl / sqrt(1 / dx^2 + 1 / dy^2) / c + + nstep = min(nstepmax, nstep) + MPI.Barrier(comm) - + # Origin of local mesh xp = coords[1] * dimx / nyp yp = coords[2] * dimy / nyp - mesh = Mesh( mx, dx, my, dy) + mesh = Mesh(mx, dx, my, dy) fields = MeshFields(mesh) - - omega = c * sqrt((md*pi/dimx)^2+(nd*pi/dimy)^2) - for j=1:my, i=1:mx - x = xp+(i-0.5)*dx - y = yp+(j-0.5)*dy - fields.bz[i,j] = (- cos(md*pi*x/dimx) - * cos(nd*pi*y/dimy) - * cos(omega*0.5*dt) ) - end - + + omega = c * sqrt((md * pi / dimx)^2 + (nd * pi / dimy)^2) + for j = 1:my, i = 1:mx + x = xp + (i - 0.5) * dx + y = yp + (j - 0.5) * dy + fields.bz[i, j] = + (-cos(md * pi * x / dimx) * cos(nd * pi * y / dimy) * cos(omega * 0.5 * dt)) + end + tag = 1111 - + for istep = 1:nstep # Loop over time - - # E(n) [1:mx]*[1:my] --> B(n+1/2) [1:mx-1]*[1:my-1] - - faraday!(fields, dt) - - - # Send to North and receive from South - MPI.Sendrecv!(view(fields.bz, 1, 1:my), north, tag, - view(fields.bz, mx+1, 1:my), south, tag, comm2d) - - # Send to West and receive from East - MPI.Sendrecv!(view(fields.bz, 1:mx, 1), west, tag, - view(fields.bz, 1:mx,my+1), east, tag, comm2d) - - # Bz(n+1/2) [1:mx]*[1:my] --> Ex(n+1) [1:mx]*[2:my+1] - # Bz(n+1/2) [1:mx]*[1:my] --> Ey(n+1) [2:mx+1]*[1:my] - - ampere_maxwell!(fields, dt) - - # Send to East and receive from West - MPI.Sendrecv!(view(fields.ex, :, my+1), east, tag, - view(fields.ex, :, 1), west, tag, comm2d) - - # Send to South and receive from North - MPI.Sendrecv!(view(fields.ey, mx+1, :), south, tag, - view(fields.ey, 1, :), north, tag, comm2d) - - # plot_fields(mesh, rank, proc, fields.bz, xp, yp, istep ) - - + + # E(n) [1:mx]*[1:my] --> B(n+1/2) [1:mx-1]*[1:my-1] + + faraday!(fields, dt) + + + # Send to North and receive from South + MPI.Sendrecv!( + view(fields.bz, 1, 1:my), + north, + tag, + view(fields.bz, mx + 1, 1:my), + south, + tag, + comm2d, + ) + + # Send to West and receive from East + MPI.Sendrecv!( + view(fields.bz, 1:mx, 1), + west, + tag, + view(fields.bz, 1:mx, my + 1), + east, + tag, + comm2d, + ) + + # Bz(n+1/2) [1:mx]*[1:my] --> Ex(n+1) [1:mx]*[2:my+1] + # Bz(n+1/2) [1:mx]*[1:my] --> Ey(n+1) [2:mx+1]*[1:my] + + ampere_maxwell!(fields, dt) + + # Send to East and receive from West + MPI.Sendrecv!( + view(fields.ex, :, my + 1), + east, + tag, + view(fields.ex, :, 1), + west, + tag, + comm2d, + ) + + # Send to South and receive from North + MPI.Sendrecv!( + view(fields.ey, mx + 1, :), + south, + tag, + view(fields.ey, 1, :), + north, + tag, + comm2d, + ) + + # plot_fields(mesh, rank, proc, fields.bz, xp, yp, istep ) + + end # next time step - + err_l2 = 0.0 - time = (nstep-0.5)*dt + time = (nstep - 0.5) * dt for j = 1:my, i = 1:mx - x = xp+(i-0.5)*dx - y = yp+(j-0.5)*dy - th_bz = (- cos(md*pi*x/dimx) - * cos(nd*pi*y/dimy) - * cos(omega*time)) - err_l2 += (fields.bz[i,j] - th_bz)^2 + x = xp + (i - 0.5) * dx + y = yp + (j - 0.5) * dy + th_bz = (-cos(md * pi * x / dimx) * cos(nd * pi * y / dimy) * cos(omega * time)) + err_l2 += (fields.bz[i, j] - th_bz)^2 end - for k in 0:proc-1 - if rank == k - println("----") - println("$rank : $mx, $my $err_l2 ") - println("$rank : $dx, $dy $err_l2 ") - end - MPI.Barrier(comm) + for k = 0:proc-1 + if rank == k + println("----") + println("$rank : $mx, $my $err_l2 ") + println("$rank : $dx, $dy $err_l2 ") + end + MPI.Barrier(comm) end - + sum_err_l2 = MPI.Allreduce(err_l2, +, comm2d) @@ -203,8 +228,8 @@ err_l2 = main(1) tend = MPI.Wtime() if MPI.Comm_rank(MPI.COMM_WORLD) == 0 - println(" error : $(err_l2) ") - println(" time : $(tend -tbegin) ") + println(" error : $(err_l2) ") + println(" time : $(tend -tbegin) ") end MPI.Finalize() diff --git a/maxwell_serial.F90 b/maxwell_serial.F90 index ac42df1..239398f 100644 --- a/maxwell_serial.F90 +++ b/maxwell_serial.F90 @@ -1,94 +1,93 @@ program Maxwell_Yee_2d -use commun -use sorties -use solveur_yee + use commun + use sorties + use solveur_yee -implicit none + implicit none -integer :: i, j -real(8) :: tcpu, x0, y0 -real(8) :: time, err_l2, th_bz -integer :: istep, iplot + integer :: i, j + real(8) :: tcpu, x0, y0 + real(8) :: time, err_l2, th_bz + integer :: istep, iplot -type(mesh_fields) :: fields, th + type(mesh_fields) :: fields, th -call cpu_time(tcpu) -call readin( ) + call cpu_time(tcpu) + call readin() -allocate(fields%ex(nx,ny+1)) -allocate(fields%ey(nx+1,ny)) -allocate(fields%bz(nx,ny)) + allocate (fields%ex(nx, ny + 1)) + allocate (fields%ey(nx + 1, ny)) + allocate (fields%bz(nx, ny)) -dx = dimx / real(nx,kind=8) -dy = dimy / real(ny,kind=8) -dt = cfl / sqrt (1./(dx*dx)+1./(dy*dy)) / c -nstep = floor(tfinal/dt) -write(*,*) -write(*,*) " dx = ", dx -write(*,*) " dy = ", dy -write(*,*) " dt = ", dt -write(*,*) -write(*,*) nx * dx -if( nstep > nstepmax ) nstep = nstepmax -write(*,*) " Nombre d'iteration nstep = ", nstep + dx = dimx/real(nx, kind=8) + dy = dimy/real(ny, kind=8) + dt = cfl/sqrt(1./(dx*dx) + 1./(dy*dy))/c + nstep = floor(tfinal/dt) + write (*, *) + write (*, *) " dx = ", dx + write (*, *) " dy = ", dy + write (*, *) " dt = ", dt + write (*, *) + write (*, *) nx*dx + if (nstep > nstepmax) nstep = nstepmax + write (*, *) " Nombre d'iteration nstep = ", nstep -time = 0. -iplot = 0 + time = 0. + iplot = 0 -omega = c * sqrt((md*pi/dimx)**2+(nd*pi/dimy)**2) + omega = c*sqrt((md*pi/dimx)**2 + (nd*pi/dimy)**2) -fields%ex(:,:) = 0.0d0 -fields%ey(:,:) = 0.0d0 + fields%ex(:, :) = 0.0d0 + fields%ey(:, :) = 0.0d0 -do j=1,ny -do i=1,nx - fields%bz(i,j) = - cos(md*pi*(i-0.5)*dx/dimx) & - * cos(nd*pi*(j-0.5)*dy/dimy) & - * cos(omega*(-0.5*dt)) -end do -end do + do j = 1, ny + do i = 1, nx + fields%bz(i, j) = -cos(md*pi*(i - 0.5)*dx/dimx) & + *cos(nd*pi*(j - 0.5)*dy/dimy) & + *cos(omega*(-0.5*dt)) + end do + end do -do istep = 1, nstep !*** Loop over time + do istep = 1, nstep !*** Loop over time - !*** Calcul de B(n+1/2) sur les pts interieurs - call faraday(fields, 1, nx, 1, ny) !Calcul de B(n-1/2)--> B(n+1/2) - - time = time + 0.5*dt + !*** Calcul de B(n+1/2) sur les pts interieurs + call faraday(fields, 1, nx, 1, ny) !Calcul de B(n-1/2)--> B(n+1/2) - call cl_periodiques(fields, 1, nx, 1, ny) + time = time + 0.5*dt - !*** Calcul de E(t=n+1) sur les pts interieurs - call ampere_maxwell(fields, 1, nx, 1, ny) + call cl_periodiques(fields, 1, nx, 1, ny) + !*** Calcul de E(t=n+1) sur les pts interieurs + call ampere_maxwell(fields, 1, nx, 1, ny) - if (mod(istep,idiag) == 0.0) then - iplot = iplot + 1 - call plot_fields(0, 1, fields, 1, nx, 1, ny, 0d0, 0d0, iplot, time ) - end if + if (mod(istep, idiag) == 0.0) then + iplot = iplot + 1 + call plot_fields(0, 1, fields, 1, nx, 1, ny, 0d0, 0d0, iplot, time) + end if - time = time + 0.5*dt + time = time + 0.5*dt -end do ! next time step + end do ! next time step -err_l2 = 0.0 -do j = 1, ny -do i = 1, nx - th_bz = - cos(md*pi*(i-0.5)*dx/dimx) & - * cos(nd*pi*(j-0.5)*dy/dimy) & - * cos(omega*(time-0.5*dt)) - err_l2 = err_l2 + (fields%bz(i,j) - th_bz)**2 -end do -end do + err_l2 = 0.0 + do j = 1, ny + do i = 1, nx + th_bz = -cos(md*pi*(i - 0.5)*dx/dimx) & + *cos(nd*pi*(j - 0.5)*dy/dimy) & + *cos(omega*(time - 0.5*dt)) + err_l2 = err_l2 + (fields%bz(i, j) - th_bz)**2 + end do + end do -call cpu_time(tcpu) + call cpu_time(tcpu) -write(*,"(10x,' istep = ',I6)",advance="no") istep -write(*,"(' time = ',e15.3,' sec')",advance="no") time -write(*,"(' erreur L2 = ',g10.5)") sqrt(err_l2) -write(*,"(//10x,' Temps CPU = ', G15.3, ' sec' )") tcpu + write (*, "(10x,' istep = ',I6)", advance="no") istep + write (*, "(' time = ',e15.3,' sec')", advance="no") time + write (*, "(' erreur L2 = ',g10.5)") sqrt(err_l2) + write (*, "(//10x,' Temps CPU = ', G15.3, ' sec' )") tcpu -stop + stop 1000 format(19f8.4) 1001 format(19i8) diff --git a/maxwell_serial.jl b/maxwell_serial.jl index c3c9c45..428a8a0 100644 --- a/maxwell_serial.jl +++ b/maxwell_serial.jl @@ -5,71 +5,75 @@ include("mesh.jl") include("fdtd.jl") -function main( nstep ) +function main(nstep) - @show cfl = 0.1 # Courant-Friedrich-Levy - @show tfinal = 10. # final time + @show cfl = 0.1 # Courant-Friedrich-Levy + @show tfinal = 10.0 # final time @show nstepmax = 1000 # max steps - @show md = 2 # md : wave number x (initial condition) - @show nd = 2 # nd : wave number y (initial condition) - @show nx = 1200 # x number of points - @show ny = 1200 # y number of points - @show dimx = 1.0 # width - @show dimy = 1.0 # height + @show md = 2 # md : wave number x (initial condition) + @show nd = 2 # nd : wave number y (initial condition) + @show nx = 1200 # x number of points + @show ny = 1200 # y number of points + @show dimx = 1.0 # width + @show dimy = 1.0 # height dx = dimx / nx dy = dimy / ny - mesh = Mesh( dimx, nx, dimy, ny) + mesh = Mesh(dimx, nx, dimy, ny) dx, dy = mesh.dx, mesh.dy - - dt = cfl / sqrt(1/dx^2+1/dy^2) / c - - @show nstep = min( nstepmax, nstep) - + + dt = cfl / sqrt(1 / dx^2 + 1 / dy^2) / c + + @show nstep = min(nstepmax, nstep) + fields = MeshFields(mesh) - - omega = c * sqrt((md*pi/dimx)^2+(nd*pi/dimy)^2) + + omega = c * sqrt((md * pi / dimx)^2 + (nd * pi / dimy)^2) # Ex and Ey are set at t = 0.0 # Bz is set at t = -dt/2 - for j=1:ny, i=1:nx - fields.bz[i,j] = (- cos(md*pi*((i-0.5)*dx/dimx)) - * cos(nd*pi*((j-0.5)*dy/dimy)) - * cos(omega*(-0.5*dt)) ) - end - + for j = 1:ny, i = 1:nx + fields.bz[i, j] = ( + -cos(md * pi * ((i - 0.5) * dx / dimx)) * + cos(nd * pi * ((j - 0.5) * dy / dimy)) * + cos(omega * (-0.5 * dt)) + ) + end + tag = 1111 - + for istep = 1:nstep # Loop over time - - # Ex(n) [1:nx]*[1:ny+1] --> B(n+1/2) [1:nx]*[1:ny] - # Ey(n) [1:nx+1]*[1:ny] --> B(n+1/2) [1:nx]*[1:ny] - - ampere_maxwell!(fields, 1, nx, 1, ny, dt) - - periodic_bc!(fields, 1, nx, 1, ny, dt) - - faraday!(fields, 1, nx, 1, ny, dt) - - err_l2 = 0.0 - time = (istep-0.5)*dt - for j = 1:ny, i = 1:nx - th_bz = (- cos(md*pi*((i-0.5)*dx/dimx)) - * cos(nd*pi*((j-0.5)*dy/dimy)) - * cos(omega*time)) - err_l2 += (fields.bz[i,j] - th_bz)^2 - end - - println(sqrt(err_l2)) + + # Ex(n) [1:nx]*[1:ny+1] --> B(n+1/2) [1:nx]*[1:ny] + # Ey(n) [1:nx+1]*[1:ny] --> B(n+1/2) [1:nx]*[1:ny] + + ampere_maxwell!(fields, 1, nx, 1, ny, dt) + + periodic_bc!(fields, 1, nx, 1, ny, dt) + + faraday!(fields, 1, nx, 1, ny, dt) + + err_l2 = 0.0 + time = (istep - 0.5) * dt + for j = 1:ny, i = 1:nx + th_bz = ( + -cos(md * pi * ((i - 0.5) * dx / dimx)) * + cos(nd * pi * ((j - 0.5) * dy / dimy)) * + cos(omega * time) + ) + err_l2 += (fields.bz[i, j] - th_bz)^2 + end + + println(sqrt(err_l2)) end # next time step - + end -main( 1 ) # trigger building +main(1) # trigger building -@time println(main( 1000 )) +@time println(main(1000)) diff --git a/mesh.jl b/mesh.jl index 052bcbf..e80c40e 100644 --- a/mesh.jl +++ b/mesh.jl @@ -1,22 +1,22 @@ struct Mesh - dimx :: Float64 - nx :: Int64 - dx :: Float64 - dimy :: Float64 - ny :: Int64 - dy :: Float64 - kx - ky + dimx::Float64 + nx::Int64 + dx::Float64 + dimy::Float64 + ny::Int64 + dy::Float64 + kx::Any + ky::Any - function Mesh( dimx, nx, dimy, ny) + function Mesh(dimx, nx, dimy, ny) dx = dimx / nx dy = dimy / ny - kx = 2π ./ dimx .* vcat(0:nx÷2-1,-nx÷2:-1) - ky = 2π ./ dimy .* vcat(0:ny÷2-1,-ny÷2:-1) |> transpose + kx = 2π ./ dimx .* vcat(0:nx÷2-1, -nx÷2:-1) + ky = 2π ./ dimy .* vcat(0:ny÷2-1, -ny÷2:-1) |> transpose - new( dimx, nx, dx, dimy, ny, dy, kx, ky) + new(dimx, nx, dx, dimy, ny, dy, kx, ky) end @@ -24,29 +24,28 @@ end struct MeshFields - mesh :: Mesh - ex :: Array{Float64, 2} - ey :: Array{Float64, 2} - bz :: Array{Float64, 2} + mesh::Mesh + ex::Array{Float64,2} + ey::Array{Float64,2} + bz::Array{Float64,2} - function MeshFields( mesh ) + function MeshFields(mesh) nx, ny = mesh.nx, mesh.ny - ex = zeros(Float64, (nx,ny+1)) - ey = zeros(Float64, (nx+1,ny)) - bz = zeros(Float64, (nx+1,ny+1)) - new( mesh, ex, ey, bz) + ex = zeros(Float64, (nx, ny + 1)) + ey = zeros(Float64, (nx + 1, ny)) + bz = zeros(Float64, (nx + 1, ny + 1)) + new(mesh, ex, ey, bz) end - function MeshFields( mesh, nx, ny ) + function MeshFields(mesh, nx, ny) - ex = zeros(Float64, (nx+1,ny+1)) - ey = zeros(Float64, (nx+1,ny+1)) - bz = zeros(Float64, (nx+1,ny+1)) - new( mesh, ex, ey, bz) + ex = zeros(Float64, (nx + 1, ny + 1)) + ey = zeros(Float64, (nx + 1, ny + 1)) + bz = zeros(Float64, (nx + 1, ny + 1)) + new(mesh, ex, ey, bz) end -end - +end diff --git a/plot_mpi.jl b/plot_mpi.jl index f7c2612..e3b67ae 100644 --- a/plot_mpi.jl +++ b/plot_mpi.jl @@ -1,49 +1,48 @@ -function plot_fields(mesh, rank, proc, field, xp, yp, iplot ) +function plot_fields(mesh, rank, proc, field, xp, yp, iplot) dx, dy = mesh.dx, mesh.dy - ix, jx = 1, mesh.nx+1 - iy, jy = 1, mesh.ny+1 + ix, jx = 1, mesh.nx + 1 + iy, jy = 1, mesh.ny + 1 if iplot == 1 mkpath("data/$rank") end io = open("data/$(rank)/$(iplot)", "w") - for j=iy:jy - for i=ix:jx - @printf( io, "%f %f %f \n", xp+(i-0.5)*dx, yp+(j-1)*dy, field[i,j]) + for j = iy:jy + for i = ix:jx + @printf(io, "%f %f %f \n", xp + (i - 0.5) * dx, yp + (j - 1) * dy, field[i, j]) end - @printf( io, "\n") + @printf(io, "\n") end close(io) - + # write master file if rank == 0 - if iplot == 1 - io = open( "bz.gnu", "w" ) - #write(io, "set xr[-0.1:1.1]\n") - #write(io, "set yr[-0.1:1.1]\n") - write(io, "set zr[-1.1:1.1]\n") - write(io, "set surf\n") - else - io = open( "bz.gnu", "a" ) - end - write(io, "set title '$(iplot)' \n") - write(io, "splot 'data/$(rank)/$(iplot)' w l") - - for p = 1:proc-1 - write(io, ", 'data/$(p)/$(iplot)' w l") - end - write( io, "\n") - #write( io, "set term gif \n") - #write( io, "set output 'image$(lpad(iplot,3,"0")).gif'\n") - #write( io, "replot\n") - - close(io) + if iplot == 1 + io = open("bz.gnu", "w") + #write(io, "set xr[-0.1:1.1]\n") + #write(io, "set yr[-0.1:1.1]\n") + write(io, "set zr[-1.1:1.1]\n") + write(io, "set surf\n") + else + io = open("bz.gnu", "a") + end + write(io, "set title '$(iplot)' \n") + write(io, "splot 'data/$(rank)/$(iplot)' w l") - end + for p = 1:proc-1 + write(io, ", 'data/$(p)/$(iplot)' w l") + end + write(io, "\n") + #write( io, "set term gif \n") + #write( io, "set output 'image$(lpad(iplot,3,"0")).gif'\n") + #write( io, "replot\n") + + close(io) -end + end +end diff --git a/read_data.jl b/read_data.jl deleted file mode 100644 index bfcfa7d..0000000 --- a/read_data.jl +++ /dev/null @@ -1,49 +0,0 @@ -open("input_data") do file - lines = readlines(file) - println(parse(Float64, lines)) -end - -IJulia.load("input_data") - -cfl = 0.2 -tfinal = 1. -nstepmax = 500 -idiag = 10 -md = 2 -nd = 2 -nx = 100 -ny = 100 -dimx = 1.0 -dimy = 1.0 - - -using JLD -jldopen("input_data.jld", "w") do file - JLD.@write file cfl - JLD.@write file tfinal - JLD.@write file nstepmax - JLD.@write file idiag - JLD.@write file md - JLD.@write file nd - JLD.@write file nx - JLD.@write file ny - JLD.@write file dimx - JLD.@write file dimy -end - - -data = load("input_data.jld") - -using JSON -open("input_data.json","w") do f - JSON.print(f, data) -end - -!cat input_data.json - -open("input_data.json", "r") do f - data = read(f, String) - println(JSON.parse(data)) -end - - diff --git a/ring.jl b/ring.jl index 554fdc1..a0b3efe 100644 --- a/ring.jl +++ b/ring.jl @@ -2,37 +2,37 @@ using MPI function main() - # Initialize the MPI environment - MPI.Init() + # Initialize the MPI environment + MPI.Init() tag = 1111 comm = MPI.COMM_WORLD - # Get the number of processes - world_size = MPI.Comm_size(comm) + # Get the number of processes + world_size = MPI.Comm_size(comm) - # Get the rank of the process - world_rank = MPI.Comm_rank(comm) + # Get the rank of the process + world_rank = MPI.Comm_rank(comm) - # Print off a hello world message - println("Processor $world_rank out of $world_size processors") + # Print off a hello world message + println("Processor $world_rank out of $world_size processors") - next_proc = mod(world_rank+1, world_size) - last_proc = mod(world_size+world_rank-1, world_size) + next_proc = mod(world_rank + 1, world_size) + last_proc = mod(world_size + world_rank - 1, world_size) - if world_rank == 0 - token = 10 # Set the token's value if you are process 0 - MPI.Send(token, next_proc, tag, comm) - (token,_) = MPI.Recv(Int, last_proc, tag, comm) - else - (token,_) = MPI.Recv(Int, last_proc, tag, comm) - MPI.Send(token+1, next_proc, tag, comm) + if world_rank == 0 + token = 10 # Set the token's value if you are process 0 + MPI.Send(token, next_proc, tag, comm) + (token, _) = MPI.Recv(Int, last_proc, tag, comm) + else + (token, _) = MPI.Recv(Int, last_proc, tag, comm) + MPI.Send(token + 1, next_proc, tag, comm) end println("Process $world_rank received token $token from process $(last_proc)") - # Finalize the MPI environment. - MPI.Finalize() + # Finalize the MPI environment. + MPI.Finalize() end diff --git a/ring_with_distributed.jl b/ring_with_distributed.jl index 0d5ed69..4cbe258 100644 --- a/ring_with_distributed.jl +++ b/ring_with_distributed.jl @@ -12,8 +12,8 @@ addprocs(4) workers() @fetch begin - println(myid()); - rand(2,2) + println(myid()) + rand(2, 2) end @sync begin @@ -26,7 +26,7 @@ println("Done!") @everywhere begin # execute this block on all workers using Random - + function complicated_calculation() sleep(1) randexp(5) # lives in Random @@ -36,14 +36,14 @@ end @fetch complicated_calculation() -ch = Channel{Int}(5) +ch = Channel{Int}(5) -isready(ch) +isready(ch) put!(ch, 3) -isready(ch) +isready(ch) take!(ch) @@ -53,10 +53,10 @@ fetch(ch) take!(ch) -isready(ch) +isready(ch) -const mychannel = RemoteChannel(()->Channel{Int}(10), workers()[2]) +const mychannel = RemoteChannel(() -> Channel{Int}(10), workers()[2]) function whohas(s::String) @everywhere begin @@ -82,7 +82,7 @@ whohas("mychannel") # + function do_something() - rc = RemoteChannel(()->Channel{Int}(10)) # lives on the master + rc = RemoteChannel(() -> Channel{Int}(10)) # lives on the master @sync for p in workers() @spawnat p put!(rc, myid()) end @@ -92,7 +92,10 @@ end r = do_something() # - -using Distributed, BenchmarkTools; rmprocs(workers()); addprocs(4); nworkers() +using Distributed, BenchmarkTools; +rmprocs(workers()); +addprocs(4); +nworkers(); # + # serial version - count heads in a series of coin tosses @@ -109,7 +112,7 @@ end # + # distributed version function add_distributed(n) - c = @distributed (+) for i in 1:n + c = @distributed (+) for i = 1:n Int(rand(Bool)) end c @@ -120,9 +123,9 @@ end # + # verbose distributed version function add_distributed(n) - c = @distributed (+) for i in 1:n + c = @distributed (+) for i = 1:n x = Int(rand(Bool)) - println(x); + println(x) x end c @@ -133,14 +136,14 @@ add_distributed(8); @everywhere using SharedArrays # must be loaded everywhere -A = rand(2,3) +A = rand(2, 3) S = SharedArray(A) # + function fill_shared_problematic(N) - S = SharedMatrix{Int64}(N,N) - @sync @distributed for i in 1:length(S) # added @sync here + S = SharedMatrix{Int64}(N, N) + @sync @distributed for i = 1:length(S) # added @sync here S[i] = i end S @@ -153,8 +156,3 @@ minimum(S) minimum(S) - - - - - diff --git a/solveur_yee.F90 b/solveur_yee.F90 index b8301c7..6eb6ea0 100644 --- a/solveur_yee.F90 +++ b/solveur_yee.F90 @@ -1,87 +1,85 @@ module solveur_yee -use commun, only: mesh_fields, dx, dy, csq , c, dt + use commun, only: mesh_fields, dx, dy, csq, c, dt -implicit none + implicit none -integer, private :: i, j -real(8), private :: dex_dx, dey_dy -real(8), private :: dex_dy, dey_dx -real(8), private :: dbz_dx, dbz_dy + integer, private :: i, j + real(8), private :: dex_dx, dey_dy + real(8), private :: dex_dy, dey_dx + real(8), private :: dbz_dx, dbz_dy contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine faraday( tm, ix, jx, iy, jy ) + subroutine faraday(tm, ix, jx, iy, jy) -type( mesh_fields ) :: tm -integer, intent(in) :: ix, jx, iy, jy + type(mesh_fields) :: tm + integer, intent(in) :: ix, jx, iy, jy !*** On utilise l'equation de Faraday sur un demi pas -!*** de temps pour le calcul du champ magnetique Bz -!*** a l'instant n puis n+1/2 +!*** de temps pour le calcul du champ magnetique Bz +!*** a l'instant n puis n+1/2 !Ex[1:nx,1:ny+1,]*Ey[1:nx+1,1:ny] -> Bz[1:nx,1:ny] -do j=iy,jy -do i=ix,jx - dex_dy = (tm%ex(i,j+1)-tm%ex(i,j)) / dy - dey_dx = (tm%ey(i+1,j)-tm%ey(i,j)) / dx - tm%bz(i,j) = tm%bz(i,j) + dt * (dex_dy - dey_dx) -end do -end do + do j = iy, jy + do i = ix, jx + dex_dy = (tm%ex(i, j + 1) - tm%ex(i, j))/dy + dey_dx = (tm%ey(i + 1, j) - tm%ey(i, j))/dx + tm%bz(i, j) = tm%bz(i, j) + dt*(dex_dy - dey_dx) + end do + end do -end subroutine faraday + end subroutine faraday !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine ampere_maxwell( tm, ix, jx, iy, jy ) + subroutine ampere_maxwell(tm, ix, jx, iy, jy) -type( mesh_fields ) :: tm -integer, intent(in) :: ix, jx, iy, jy + type(mesh_fields) :: tm + integer, intent(in) :: ix, jx, iy, jy !*** Calcul du champ electrique E au temps n+1 !*** sur les points internes du maillage !*** Ex aux points (i+1/2,j) !*** Ey aux points (i,j+1/2) -do j=iy+1,jy -do i=ix,jx - dbz_dy = (tm%bz(i,j)-tm%bz(i,j-1)) / dy - tm%ex(i,j) = tm%ex(i,j) + dt*csq*dbz_dy -end do -end do + do j = iy + 1, jy + do i = ix, jx + dbz_dy = (tm%bz(i, j) - tm%bz(i, j - 1))/dy + tm%ex(i, j) = tm%ex(i, j) + dt*csq*dbz_dy + end do + end do -do j=iy,jy -do i=ix+1,jx - dbz_dx = (tm%bz(i,j)-tm%bz(i-1,j)) / dx - tm%ey(i,j) = tm%ey(i,j) - dt*csq*dbz_dx -end do -end do + do j = iy, jy + do i = ix + 1, jx + dbz_dx = (tm%bz(i, j) - tm%bz(i - 1, j))/dx + tm%ey(i, j) = tm%ey(i, j) - dt*csq*dbz_dx + end do + end do -end subroutine ampere_maxwell + end subroutine ampere_maxwell !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine cl_periodiques(tm, ix, jx, iy, jy) + subroutine cl_periodiques(tm, ix, jx, iy, jy) -type( mesh_fields ) :: tm -integer, intent(in) :: ix, jx, iy, jy + type(mesh_fields) :: tm + integer, intent(in) :: ix, jx, iy, jy + do i = ix, jx + dbz_dy = (tm%bz(i, iy) - tm%bz(i, jy))/dy + tm%ex(i, iy) = tm%ex(i, iy) + dt*csq*dbz_dy + tm%ex(i, jy + 1) = tm%ex(i, iy) + end do -do i = ix, jx - dbz_dy = (tm%bz(i,iy)-tm%bz(i,jy)) / dy - tm%ex(i,iy) = tm%ex(i,iy) + dt*csq*dbz_dy - tm%ex(i,jy+1) = tm%ex(i,iy) -end do + do j = iy, jy + dbz_dx = (tm%bz(ix, j) - tm%bz(jx, j))/dx + tm%ey(ix, j) = tm%ey(ix, j) - dt*csq*dbz_dx + tm%ey(jx + 1, j) = tm%ey(ix, j) + end do - -do j = iy, jy - dbz_dx = (tm%bz(ix,j)-tm%bz(jx,j)) / dx - tm%ey(ix,j) = tm%ey(ix,j) - dt*csq*dbz_dx - tm%ey(jx+1,j) = tm%ey(ix,j) -end do - -end subroutine + end subroutine end module solveur_yee diff --git a/sorties.F90 b/sorties.F90 index 92dbacb..d6048ea 100644 --- a/sorties.F90 +++ b/sorties.F90 @@ -1,80 +1,80 @@ module sorties -use commun, only: mesh_fields, dx, dy + use commun, only: mesh_fields, dx, dy -implicit none + implicit none -integer, private :: i, j, k -character(len=10) :: outdir = 'data/' + integer, private :: i, j, k + character(len=10) :: outdir = 'data/' contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine plot_fields(rang, nproc, f, ix, jx, iy, jy, & - xp, yp, iplot, time ) + subroutine plot_fields(rang, nproc, f, ix, jx, iy, jy, & + xp, yp, iplot, time) -integer, intent(in) :: rang, ix, jx, iy, jy, nproc -type (mesh_fields), intent(in) :: f -integer :: iplot, i, j -real(8) :: time, xp, yp -integer :: kk0, kk1, kk2, kk3, kk4 -character(len=4) :: fin -character(len=1) :: aa,bb,cc,dd + integer, intent(in) :: rang, ix, jx, iy, jy, nproc + type(mesh_fields), intent(in) :: f + integer :: iplot, i, j + real(8) :: time, xp, yp + integer :: kk0, kk1, kk2, kk3, kk4 + character(len=4) :: fin + character(len=1) :: aa, bb, cc, dd -if (iplot == 1) then - k = len_trim(outdir) - outdir(k:) = "/"//char(rang+48)//"/" - call system("mkdir -p "//outdir) -end if + if (iplot == 1) then + k = len_trim(outdir) + outdir(k:) = "/"//char(rang + 48)//"/" + call system("mkdir -p "//outdir) + end if -kk0 = iplot -kk1 = kk0/1000 -aa = char(kk1 + 48) -kk2 = (kk0 - kk1*1000)/100 -bb = char(kk2 + 48) -kk3 = (kk0 - (kk1*1000) - (kk2*100))/10 -cc = char(kk3 + 48) -kk4 = (kk0 - (kk1*1000) - (kk2*100) - (kk3*10))/1 -dd = char(kk4 + 48) -fin = aa//bb//cc//dd + kk0 = iplot + kk1 = kk0/1000 + aa = char(kk1 + 48) + kk2 = (kk0 - kk1*1000)/100 + bb = char(kk2 + 48) + kk3 = (kk0 - (kk1*1000) - (kk2*100))/10 + cc = char(kk3 + 48) + kk4 = (kk0 - (kk1*1000) - (kk2*100) - (kk3*10))/1 + dd = char(kk4 + 48) + fin = aa//bb//cc//dd -open( 80, file = trim(outdir)//fin//".dat" ) - do i=ix,jx - do j=iy,jy - write(80,"(4e10.2)") xp+(i-0.5)*dx, yp+(j-.5)*dy, f%bz(i,j) + open (80, file=trim(outdir)//fin//".dat") + do i = ix, jx + do j = iy, jy + write (80, "(4e10.2)") xp + (i - 0.5)*dx, yp + (j - .5)*dy, f%bz(i, j) + end do + write (80, *) end do - write(80,*) - end do -close(80) - -if (rang == 0) then - do k = 1, 3 - open( 90, file = 'bz.gnu', position="append" ) - if ( iplot == 1 ) then - rewind(90) - write(90,*)"set xr[-0.1:1.1]" - write(90,*)"set yr[-0.1:1.1]" - write(90,*)"set zr[-1.1:1.1]" - !write(90,*)"set cbrange[-1:1]" - !write(90,*)"set pm3d" - write(90,*)"set surf" - !write(90,*)"set term x11" + close (80) + + if (rang == 0) then + do k = 1, 3 + open (90, file='bz.gnu', position="append") + if (iplot == 1) then + rewind (90) + write (90, *) "set xr[-0.1:1.1]" + write (90, *) "set yr[-0.1:1.1]" + write (90, *) "set zr[-1.1:1.1]" + !write(90,*)"set cbrange[-1:1]" + !write(90,*)"set pm3d" + write (90, *) "set surf" + !write(90,*)"set term x11" + end if + write (90, *) "set title 'Time = ", time, "'" + write (90, "(a)", advance='no') "splot '" & + & //trim(outdir)//fin//".dat' u 1:2:3 w lines" + + do j = 1, nproc - 1 + write (90, "(a)", advance='no') "& + & ,'"//outdir(1:5)//char(j + 48)//"/"//fin// & + & ".dat' u 1:2:3 w lines" + end do + write (90, *) + close (90) + end do end if - write(90,*)"set title 'Time = ",time,"'" - write(90,"(a)",advance='no')"splot '" & - & //trim(outdir)//fin//".dat' u 1:2:3 w lines" - - do j = 1, nproc - 1 - write(90,"(a)",advance='no')"& - & ,'"//outdir(1:5)//char(j+48)//"/"//fin// & - & ".dat' u 1:2:3 w lines" - end do - write(90,*) - close(90) - end do -end if -end subroutine plot_fields + end subroutine plot_fields end module sorties