Skip to content

Commit f94252a

Browse files
loiseaujcloiseaujc
authored
Modernized solver (#8)
Merci !
1 parent 9120655 commit f94252a

File tree

5 files changed

+322
-0
lines changed

5 files changed

+322
-0
lines changed

‎modern_fortran/Makefile‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
SRCS = solveur_yee.f08 sorties.f08 commun.f08
2+
OBJS = solveur_yee.o sorties.o commun.o
3+
F90 = mpif90
4+
OPT = -O3 -mtune=native -march=native
5+
F90FLAGS = -cpp $(OPT)
6+
LDFLAGS = $(OPT)
7+
8+
all: maxwell_serial_fortran
9+
10+
maxwell_serial_fortran: maxwell_serial.o
11+
$(F90) $(LDFLAGS) -o $@ $(OBJS) $^ $(LIBS)
12+
13+
clean:
14+
rm -rf $(PROG) *.o *.mod data/* *.gnu core error.* output.* *.a
15+
16+
.SUFFIXES: $(SUFFIXES) .f08
17+
18+
.f08.o:
19+
$(F90) $(F90FLAGS) -c $<
20+
21+
.mod.o:
22+
$(F90) $(F90FLAGS) -c $*.f08
23+
24+
commun.o: commun.f08
25+
solveur_yee.o: solveur_yee.f08 commun.o
26+
sorties.o: sorties.f08 commun.o
27+
maxwell_serial.o: maxwell_serial.f08 $(OBJS)

‎modern_fortran/commun.f08‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
module commun
2+
use, intrinsic:: iso_fortran_env, only: dp => real64, ilp => int32
3+
implicit none(external)
4+
public
5+
6+
type mesh_fields
7+
real(dp), dimension(:, :), pointer :: ex, ey
8+
real(dp), dimension(:, :), pointer :: bz
9+
end type mesh_fields
10+
11+
real(dp), parameter :: c = 1.0_dp, csq = c**2
12+
real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp)
13+
14+
integer(ilp) :: md, nd
15+
integer(ilp) :: nx, ny
16+
integer(ilp) :: nstep, nstepmax
17+
integer(ilp) :: idiag
18+
19+
real(dp) :: dt, dx, dy
20+
real(dp) :: dimx, dimy
21+
real(dp) :: cfl
22+
real(dp) :: tfinal
23+
real(dp) :: omega
24+
25+
contains
26+
27+
subroutine readin()
28+
implicit none(external)
29+
namelist /donnees/ cfl, & !nbre de Courant
30+
tfinal, & !duree maxi
31+
nstepmax, & !nbre d'iterations maxi
32+
idiag, & !frequence des diagnostics
33+
md, nd, & !nombre d'onde de la solution initiale
34+
nx, ny, & !nombre de points dans les deux directions
35+
dimx, dimy !dimensions du domaine
36+
37+
!***Initialisation des valeurs pas default
38+
open (10, file="../input_data", status='old')
39+
read (10, donnees)
40+
close (10)
41+
42+
write (*, "(a,g12.3)") " largeur dimx = ", dimx
43+
write (*, "(a,g12.3)") " longueur dimy = ", dimy
44+
write (*, "(a,g12.3)") " temps = ", tfinal
45+
write (*, "(a,g12.3)") " nombre nx = ", nx
46+
write (*, "(a,g12.3)") " nombre ny = ", ny
47+
write (*, "(a,g12.3)") " nombre de Courant = ", cfl
48+
write (*, "(a,g12.3)") " frequence des diagnostics = ", idiag
49+
end subroutine readin
50+
51+
end module commun

‎modern_fortran/maxwell_serial.f08‎

Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
program Maxwell_Yee_2d
2+
use, intrinsic :: iso_fortran_env, only: output_unit
3+
use commun
4+
use sorties, only: plot_fields
5+
use solveur_yee, only: faraday, ampere_maxwell, cl_periodiques
6+
implicit none(external)
7+
8+
integer(ilp) :: i, j
9+
real(dp) :: tcpu, x0, y0
10+
real(dp) :: time, err_l2, th_bz
11+
integer(ilp) :: istep, iplot
12+
type(mesh_fields) :: fields, th
13+
14+
call cpu_time(tcpu)
15+
call readin()
16+
17+
allocate (fields%ex(nx, ny + 1), source=0.0_dp)
18+
allocate (fields%ey(nx + 1, ny), source=0.0_dp)
19+
allocate (fields%bz(nx, ny))
20+
21+
dx = dimx/real(nx, kind=dp)
22+
dy = dimy/real(ny, kind=dp)
23+
dt = cfl/sqrt(1.0_dp/(dx*dx) + 1.0_dp/(dy*dy))/c
24+
nstep = floor(tfinal/dt)
25+
write (output_unit, *)
26+
write (output_unit, *) " dx = ", dx
27+
write (output_unit, *) " dy = ", dy
28+
write (output_unit, *) " dt = ", dt
29+
write (output_unit, *)
30+
write (output_unit, *) nx*dx
31+
nstep = min(nstep, nstepmax)
32+
write (output_unit, *) " Nombre d'iteration nstep = ", nstep
33+
34+
time = 0.0_dp
35+
iplot = 0
36+
37+
omega = c*sqrt((md*pi/dimx)**2 + (nd*pi/dimy)**2)
38+
39+
do concurrent(i=1:nx, j=1:ny)
40+
fields%bz(i, j) = -cos(md*pi*(i - 0.5_dp)*dx/dimx) &
41+
*cos(nd*pi*(j - 0.5_dp)*dy/dimy) &
42+
*cos(omega*(-0.5_dp*dt))
43+
end do
44+
45+
do istep = 1, nstep !*** Loop over time
46+
47+
!*** Calcul de B(n+1/2) sur les pts interieurs
48+
call faraday(fields, 1, nx, 1, ny) !Calcul de B(n-1/2)--> B(n+1/2)
49+
50+
time = time + 0.5_dp*dt
51+
52+
call cl_periodiques(fields, 1, nx, 1, ny)
53+
54+
!*** Calcul de E(t=n+1) sur les pts interieurs
55+
call ampere_maxwell(fields, 1, nx, 1, ny)
56+
57+
if (mod(istep, idiag) == 0) then
58+
iplot = iplot + 1
59+
call plot_fields(0, 1, fields, 1, nx, 1, ny, 0.0_dp, 0.0_dp, iplot, time)
60+
end if
61+
62+
time = time + 0.5_dp*dt
63+
64+
end do ! next time step
65+
66+
err_l2 = 0.0_dp
67+
do concurrent(i=1:nx, j=1:ny) reduce(+:err_l2) local(th_bz)
68+
th_bz = -cos(md*pi*(i - 0.5)*dx/dimx) &
69+
*cos(nd*pi*(j - 0.5)*dy/dimy) &
70+
*cos(omega*(time - 0.5*dt))
71+
err_l2 = err_l2 + (fields%bz(i, j) - th_bz)**2
72+
end do
73+
74+
call cpu_time(tcpu)
75+
76+
write (output_unit, "(10x,' istep = ',I6)", advance="no") istep
77+
write (output_unit, "(' time = ',e15.3,' sec')", advance="no") time
78+
write (output_unit, "(' erreur L2 = ',g10.5)") sqrt(err_l2)
79+
write (output_unit, "(//10x,' Temps CPU = ', G15.3, ' sec' )") tcpu
80+
81+
1000 format(19f8.4)
82+
1001 format(19i8)
83+
84+
end program Maxwell_Yee_2d

‎modern_fortran/solveur_yee.f08‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
module solveur_yee
2+
use commun, only: dp, ilp, mesh_fields, dx, dy, csq, c, dt
3+
implicit none(external)
4+
private
5+
6+
interface
7+
pure module subroutine faraday(tm, ix, jx, iy, jy)
8+
implicit none(external)
9+
type(mesh_fields), intent(inout) :: tm
10+
integer(ilp), intent(in) :: ix, jx, iy, jy
11+
end subroutine faraday
12+
pure module subroutine ampere_maxwell(tm, ix, jx, iy, jy)
13+
implicit none(external)
14+
type(mesh_fields), intent(inout) :: tm
15+
integer(ilp), intent(in) :: ix, jx, iy, jy
16+
end subroutine ampere_maxwell
17+
pure module subroutine cl_periodiques(tm, ix, jx, iy, jy)
18+
implicit none(external)
19+
type(mesh_fields), intent(inout) :: tm
20+
integer(ilp), intent(in) :: ix, jx, iy, jy
21+
end subroutine cl_periodiques
22+
end interface
23+
public :: faraday, ampere_maxwell, cl_periodiques
24+
25+
contains
26+
27+
module procedure faraday
28+
integer(ilp) :: i, j
29+
real(dp) :: dex_dy, dey_dx
30+
!*** On utilise l'equation de Faraday sur un demi pas
31+
!*** de temps pour le calcul du champ magnetique Bz
32+
!*** a l'instant n puis n+1/2
33+
!Ex[1:nx,1:ny+1,]*Ey[1:nx+1,1:ny] -> Bz[1:nx,1:ny]
34+
do concurrent(i=ix:jx, j=iy:jy) local(dex_dy, dey_dx)
35+
dex_dy = (tm%ex(i, j + 1) - tm%ex(i, j))/dy
36+
dey_dx = (tm%ey(i + 1, j) - tm%ey(i, j))/dx
37+
tm%bz(i, j) = tm%bz(i, j) + dt*(dex_dy - dey_dx)
38+
end do
39+
end procedure faraday
40+
41+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42+
43+
module procedure ampere_maxwell
44+
integer(ilp) :: i, j
45+
real(dp) :: dbz_dy, dbz_dx
46+
!*** Calcul du champ electrique E au temps n+1
47+
!*** sur les points internes du maillage
48+
!*** Ex aux points (i+1/2,j)
49+
!*** Ey aux points (i,j+1/2)
50+
do concurrent(i=ix:jx, j=iy + 1:jy) local(dbz_dy)
51+
dbz_dy = (tm%bz(i, j) - tm%bz(i, j - 1))/dy
52+
tm%ex(i, j) = tm%ex(i, j) + dt*csq*dbz_dy
53+
end do
54+
do concurrent(i=ix + 1:jx, j=iy:jy) local(dbz_dx)
55+
dbz_dx = (tm%bz(i, j) - tm%bz(i - 1, j))/dx
56+
tm%ey(i, j) = tm%ey(i, j) - dt*csq*dbz_dx
57+
end do
58+
end procedure ampere_maxwell
59+
60+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
61+
62+
module procedure cl_periodiques
63+
integer(ilp) :: i, j
64+
real(dp) :: dbz_dy, dbz_dx
65+
do concurrent(i=ix:jx) local(dbz_dy)
66+
dbz_dy = (tm%bz(i, iy) - tm%bz(i, jy))/dy
67+
tm%ex(i, iy) = tm%ex(i, iy) + dt*csq*dbz_dy
68+
tm%ex(i, jy + 1) = tm%ex(i, iy)
69+
end do
70+
do concurrent(j=iy:jy) local(dbz_dx)
71+
dbz_dx = (tm%bz(ix, j) - tm%bz(jx, j))/dx
72+
tm%ey(ix, j) = tm%ey(ix, j) - dt*csq*dbz_dx
73+
tm%ey(jx + 1, j) = tm%ey(ix, j)
74+
end do
75+
end procedure cl_periodiques
76+
77+
end module solveur_yee

‎modern_fortran/sorties.f08‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
module sorties
2+
use commun, only: dp, ilp, mesh_fields, dx, dy
3+
implicit none(external)
4+
private
5+
6+
character(len=10) :: outdir = '../data/'
7+
8+
interface
9+
module subroutine plot_fields(rang, nproc, f, ix, jx, iy, jy, &
10+
xp, yp, iplot, time)
11+
implicit none(external)
12+
integer(ilp), intent(in) :: rang, ix, jx, iy, jy, nproc, iplot
13+
type(mesh_fields), intent(in) :: f
14+
real(dp), intent(in) :: xp, yp, time
15+
end subroutine plot_fields
16+
end interface
17+
public :: plot_fields
18+
19+
contains
20+
21+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
22+
module procedure plot_fields
23+
integer :: i, j, k
24+
integer :: kk0, kk1, kk2, kk3, kk4
25+
character(len=4) :: fin
26+
character(len=1) :: aa, bb, cc, dd
27+
28+
if (iplot == 1) then
29+
k = len_trim(outdir)
30+
outdir(k:) = "/"//char(rang + 48)//"/"
31+
call system("mkdir -p "//outdir)
32+
end if
33+
34+
kk0 = iplot
35+
kk1 = kk0/1000
36+
aa = char(kk1 + 48)
37+
kk2 = (kk0 - kk1*1000)/100
38+
bb = char(kk2 + 48)
39+
kk3 = (kk0 - (kk1*1000) - (kk2*100))/10
40+
cc = char(kk3 + 48)
41+
kk4 = (kk0 - (kk1*1000) - (kk2*100) - (kk3*10))/1
42+
dd = char(kk4 + 48)
43+
fin = aa//bb//cc//dd
44+
45+
open (80, file=trim(outdir)//fin//".dat")
46+
do i = ix, jx
47+
do j = iy, jy
48+
write (80, "(4e10.2)") xp + (i - 0.5)*dx, yp + (j - .5)*dy, f%bz(i, j)
49+
end do
50+
write (80, *)
51+
end do
52+
close (80)
53+
54+
if (rang == 0) then
55+
do k = 1, 3
56+
open (90, file='bz.gnu', position="append")
57+
if (iplot == 1) then
58+
rewind (90)
59+
write (90, *) "set xr[-0.1:1.1]"
60+
write (90, *) "set yr[-0.1:1.1]"
61+
write (90, *) "set zr[-1.1:1.1]"
62+
!write(90,*)"set cbrange[-1:1]"
63+
!write(90,*)"set pm3d"
64+
write (90, *) "set surf"
65+
!write(90,*)"set term x11"
66+
end if
67+
write (90, *) "set title 'Time = ", time, "'"
68+
write (90, "(a)", advance='no') "splot '" &
69+
& //trim(outdir)//fin//".dat' u 1:2:3 w lines"
70+
71+
do j = 1, nproc - 1
72+
write (90, "(a)", advance='no') "&
73+
& ,'"//outdir(1:5)//char(j + 48)//"/"//fin// &
74+
& ".dat' u 1:2:3 w lines"
75+
end do
76+
write (90, *)
77+
close (90)
78+
end do
79+
end if
80+
81+
end procedure plot_fields
82+
83+
end module sorties

0 commit comments

Comments
 (0)