From 20f590fda2a58715e8fdbb9422a8f1b140cf6490 Mon Sep 17 00:00:00 2001 From: Pierre Navaro Date: Tue, 19 Mar 2019 16:16:05 +0100 Subject: [PATCH] Implement initialization of Maxwell solver 1d with FEM --- .gitignore | 2 + Project.toml | 9 + src/GEMPIC.jl | 5 + src/maxwell_1d_fem.jl | 419 ++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 177 ++++++++++++++++++ 5 files changed, 612 insertions(+) create mode 100644 Project.toml create mode 100644 src/GEMPIC.jl create mode 100644 src/maxwell_1d_fem.jl create mode 100644 test/runtests.jl diff --git a/.gitignore b/.gitignore index 381e0b6..e09d7db 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ *.jl.*.cov *.jl.mem deps/deps.jl +Manifest.toml +*.swp diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..a400641 --- /dev/null +++ b/Project.toml @@ -0,0 +1,9 @@ +name = "GEMPIC" +uuid = "b6d65c3a-4a4e-11e9-25d0-d309dc85ddeb" +authors = ["Pierre Navaro "] +version = "0.1.0" + +[deps] +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/GEMPIC.jl b/src/GEMPIC.jl new file mode 100644 index 0000000..d474359 --- /dev/null +++ b/src/GEMPIC.jl @@ -0,0 +1,5 @@ +module GEMPIC + +include("maxwell_1d_fem.jl") + +end diff --git a/src/maxwell_1d_fem.jl b/src/maxwell_1d_fem.jl new file mode 100644 index 0000000..dea1a8c --- /dev/null +++ b/src/maxwell_1d_fem.jl @@ -0,0 +1,419 @@ +using FFTW + +export Maxwell1DFEM + +""" + maxwell_solver = MaxwellFEM1D( domain, ncells, degree ) + +1D Maxwell spline finite element solver on a periodic grid + +- Lx : length of Periodic domain +- delta_x : cell size +- n_dofs : number of cells (and grid points) +- s_deg_0 : spline degree 0-forms +- s_deg_1 : spline degree 1-forms +- mass_0 : coefficients of 0-form mass matrix +- mass_1 : coefficients of 1-form mass matrix +- eig_mass0 : eigenvalues of circulant 0-form mass matrix +- eig_mass1 : eigenvalues of circulant 1-form mass matrix +- eig_weak_ampere : eigenvalues of circulant update matrix for Ampere +- eig_weak_poisson : eigenvalues of circulant update matrix for Poisson +- plan_fw : fft plan (forward) +- plan_bw : fft plan (backward) + +""" +mutable struct Maxwell1DFEM + + Lx :: Float64 + delta_x :: Float64 + n_dofs :: Int32 + s_deg_0 :: Int32 + s_deg_1 :: Int32 + mass_0 :: Vector{Float64} + mass_1 :: Vector{Float64} + eig_mass0 :: Vector{Float64} + eig_mass1 :: Vector{Float64} + eig_weak_ampere :: Vector{Float64} + eig_weak_poisson :: Vector{Float64} + + plan_fw :: FFTW.FFTWPlan + plan_bw :: FFTW.FFTWPlan + work :: Vector{Float64} + wsave :: Vector{Float64} + + function Maxwell1DFEM( domain, ncells :: Int, degree :: Int ) + + n_dofs = ncells + Lx = domain[2] - domain[1] + delta_x = Lx / n_dofs + s_deg_0 = degree + s_deg_1 = degree - 1 + + mass_0 = zeros(Float64, s_deg_0+1) + mass_1 = zeros(Float64, s_deg_0) + + if s_deg_0 == 1 # linear and constant splines + # Upper diagonal coeficients of linear spline mass matrix (from Eulerian numbers) + mass_0[1] = 4.0/6.0 + mass_0[2] = 1.0/6.0 + # Upper diagonal coeficients of constant spline mass matrix + mass_1[1] = 1.0 + elseif s_deg_0 == 2 # quadratic and linear splines + # Upper diagonal coeficients of quadratic spline mass matrix (from Eulerian numbers) + mass_0[1] = 66.0/120.0 + mass_0[2] = 26.0/120.0 + mass_0[3] = 1.0/120.0 + # Upper diagonal coeficients of linear spline mass matrix (from Eulerian numbers) + mass_1[1] = 4.0/6.0 + mass_1[2] = 1.0/6.0 + elseif s_deg_0 == 3 + # Upper diagonal coeficients of cubic spline mass matrix (from Eulerian numbers) + mass_0[1] = 2416.0/5040.0 + mass_0[2] = 1191.0/5040.0 + mass_0[3] = 120.0/5040.0 + mass_0[4] = 1.0/5040.0 + # Upper diagonal coeficients of quadratic spline mass matrix (from Eulerian numbers) + mass_1[1] = 66.0/120.0 + mass_1[2] = 26.0/120.0 + mass_1[3] = 1.0/120.0 + else + throw( ArgumentError("Wrong value of degree = $degree") ) + end + + eig_mass0 = zeros(Float64, n_dofs) + eig_mass1 = zeros(Float64, n_dofs) + eig_weak_ampere = zeros(Float64, n_dofs) + eig_weak_poisson = zeros(Float64, n_dofs) + + work = zeros(Float64, n_dofs) + wsave = zeros(Float64, n_dofs) + plan_fw = FFTW.plan_r2r(work, FFTW.R2HC) + plan_bw = FFTW.plan_r2r(wsave, FFTW.HC2R) + + # Compute eigenvalues of circulant Ampere update matrix M_0^{-1} D^T M_1 + # and circulant Poisson Matrix (D^T M_1 D)^{-1} + # zero mode vanishes due to derivative matrix D^T + + eig_weak_ampere[1] = 0.0 + eig_weak_poisson[1] = 0.0 # Matrix is not invertible: 0-mode is set to 0 + eig_mass0[1] = 1.0 # sum of coefficents is one + eig_mass1[1] = 1.0 # sum of coefficents is one + + for k = 1:n_dofs÷2-1 + + coef0 = mass_0[1] + coef1 = mass_1[1] + for j = 1:s_deg_0-1 + cos_mode = cos(2*pi*j*k÷n_dofs) + coef0 = coef0 + 2 * mass_0[j+1] * cos_mode + coef1 = coef1 + 2 * mass_1[j+1] * cos_mode + end + # add last term for larger matrix + j = s_deg_0 + coef0 = coef0 + 2 * mass_0[j+1]*cos(2*pi*j*k÷n_dofs) + # compute eigenvalues + eig_mass0[k+1] = coef0 # real part + eig_mass0[n_dofs-k+1] = 0.0 # imaginary part + eig_mass1[k+1] = coef1 # real part + eig_mass1[n_dofs-k+1] = 0.0 # imaginary part + cos_mode = cos(2*pi*k÷n_dofs) + sin_mode = sin(2*pi*k÷n_dofs) + eig_weak_ampere[k+1] = (coef1 / coef0) * (1-cos_mode) # real part + eig_weak_ampere[n_dofs-k+1] = -(coef1 / coef0) * sin_mode # imaginary part + eig_weak_poisson[k+1] = 1.0 / (coef1 * ((1-cos_mode)^2 + sin_mode^2)) # real part + eig_weak_poisson[n_dofs-k+1] = 0.0 # imaginary part + + end + + # N/2 mode + coef0 = mass_0[1] + coef1 = mass_1[1] + for j=1:s_deg_0 - 1 + coef0 = coef0 + 2 * mass_0[j+1]*cos(pi*j) + coef1 = coef1 + 2 * mass_1[j+1]*cos(pi*j) + end + + # add last term for larger matrix + j = s_deg_0 + coef0 = coef0 + 2 * mass_0[j+1]*cos(pi*j) + + # compute eigenvalues + eig_mass0[n_dofs÷2+1] = coef0 + eig_mass1[n_dofs÷2+1] = coef1 + eig_weak_ampere[n_dofs÷2+1] = 2.0 * (coef1 / coef0) + eig_weak_poisson[n_dofs÷2+1] = 1.0 / (coef1 * 4.0) + + new( Lx, delta_x, n_dofs, s_deg_0, s_deg_1, mass_0, mass_1, + eig_mass0, eig_mass1, eig_weak_ampere, eig_weak_poisson, + plan_fw, plan_bw, work, wsave ) + + + end + +end + +#= + +contains + !> compute Ey from Bz using weak Ampere formulation + subroutine sll_s_compute_e_from_b_1d_fem(self, delta_t, field_in, field_out) + class(sll_t_maxwell_1d_fem) :: self + sll_real64, intent(in) :: delta_t !< Time step + sll_real64, intent(in) :: field_in(:) !< Bz + sll_real64, intent(inout) :: field_out(:) !< Ey + ! local variables + sll_real64 :: coef + + ! Compute potential weak curl of Bz using eigenvalue of circulant inverse matrix + call solve_circulant(self, self%eig_weak_ampere, field_in, self%work) + ! Update bz from self value + coef = delta_t/self%delta_x + field_out = field_out + coef*self%work + + end subroutine sll_s_compute_e_from_b_1d_fem + + !> Compute Bz from Ey using strong 1D Faraday equation for spline coefficients + !> $B_z^{new}(x_j) = B_z^{old}(x_j) - \frac{\Delta t}{\Delta x} (E_y(x_j) - E_y(x_{j-1}) $ + subroutine sll_s_compute_b_from_e_1d_fem(self, delta_t, field_in, field_out) + class(sll_t_maxwell_1d_fem) :: self + sll_real64, intent(in) :: delta_t + sll_real64, intent(in) :: field_in(:) ! ey + sll_real64, intent(inout) :: field_out(:) ! bz + ! local variables + sll_real64 :: coef + sll_int32 :: i + + coef = delta_t/self%delta_x + ! relation betwen spline coefficients for strong Ampere + do i=2,self%n_dofs + field_out(i) = field_out(i) + coef * ( field_in(i-1) - field_in(i) ) + end do + ! treat Periodic point + field_out(1) = field_out(1) + coef * ( field_in(self%n_dofs) - field_in(1) ) + end subroutine sll_s_compute_b_from_e_1d_fem + + !> Compute E_i from j_i integrated over the time interval using weak Ampere formulation + subroutine compute_E_from_j_1d_fem(self, current, component, E) + class(sll_t_maxwell_1d_fem) :: self !< Maxwell solver class + sll_real64,dimension(:),intent(in) :: current !< Component \a component of the current integrated over time interval + sll_int32, intent(in) :: component !< Component of the Efield to be computed + sll_real64,dimension(:),intent(inout) :: E !< Updated electric field + ! local variables + sll_int32 :: i + sll_real64, dimension(self%n_dofs) :: eigvals + + ! Multiply by inverse mass matrix using the eigenvalues of the circulant inverse matrix + eigvals=0.0 + if (component == 1) then + do i=1,self%n_dofs/2+1 + eigvals(i) = 1.0 / self%eig_mass1(i) + end do + call solve_circulant(self, eigvals, current, self%work) + elseif (component == 2) then + do i=1,self%n_dofs/2+1 + eigvals(i) = 1.0 / self%eig_mass0(i) + end do + call solve_circulant(self, eigvals, current, self%work) + else + print*, 'Component ', component, 'not implemented in compute_E_from_j_1d_fem.' + end if + + + ! Update the electric field and scale + E = E - self%work/self%delta_x + + end subroutine compute_E_from_j_1d_fem + + subroutine sll_s_compute_e_from_rho_1d_fem(self, E, rho ) + class(sll_t_maxwell_1d_fem) :: self + sll_real64,dimension(:),intent(in) :: rho + sll_real64,dimension(:),intent(out) :: E + ! local variables + sll_int32 :: i + + ! Compute potential phi from rho, using eigenvalue of circulant inverse matrix + call solve_circulant(self, self%eig_weak_poisson, rho, self%work) + ! Compute spline coefficients of Ex from those of phi + do i=2,self%n_dofs + E(i) = (self%work(i-1) - self%work(i)) !* (self%delta_x) + end do + ! treat Periodic point + E(1) = (self%work(self%n_dofs) - self%work(1)) !* (self%delta_x) + + end subroutine sll_s_compute_e_from_rho_1d_fem + + subroutine solve_circulant(self, eigvals, rhs, res) + class(sll_t_maxwell_1d_fem) :: self + sll_real64, intent(in) :: eigvals(:) ! eigenvalues of circulant matrix + sll_real64, intent(in) :: rhs(:) + sll_real64, intent(out) :: res(:) + ! local variables + sll_int32 :: k + sll_real64 :: re, im + + ! Compute res from rhs, using eigenvalue of circulant matrix + res = rhs + ! Forward FFT + call sll_s_fft_exec_r2r_1d ( self%plan_fw, res, self%wsave ) + self%wsave(1) = self%wsave(1) * eigvals(1) + do k=2, self%n_dofs/2 + re = self%wsave(k) * eigvals(k) - & + self%wsave(self%n_dofs-k+2) * eigvals(self%n_dofs-k+2) + im = self%wsave(k) * eigvals(self%n_dofs-k+2) + & + self%wsave(self%n_dofs-k+2) * eigvals(k) + self%wsave(k) = re + self%wsave(self%n_dofs-k+2) = im + end do + self%wsave(self%n_dofs/2+1) = self%wsave(self%n_dofs/2+1)*eigvals(self%n_dofs/2+1) + ! Backward FFT + call sll_s_fft_exec_r2r_1d( self%plan_bw, self%wsave, res ) + ! normalize + res = res / self%n_dofs + end subroutine solve_circulant + + + !> Compute the FEM right-hand-side for a given function f and periodic splines of given degree + !> Its components are $\int f N_i dx$ where $N_i$ is the B-spline starting at $x_i$ + subroutine sll_s_compute_fem_rhs(self, func, degree, coefs_dofs) + class(sll_t_maxwell_1d_fem) :: self + procedure(sll_i_function_1d_real64) :: func + sll_int32, intent(in) :: degree + sll_real64, intent(out) :: coefs_dofs(:) ! Finite Element right-hand-side + ! local variables + sll_int32 :: i,j,k + sll_real64 :: coef + sll_real64, dimension(2,degree+1) :: xw_gauss + sll_real64, dimension(degree+1,degree+1) :: bspl + + ! take enough Gauss points so that projection is exact for splines of degree deg + ! rescale on [0,1] for compatibility with B-splines + xw_gauss = sll_f_gauss_legendre_points_and_weights(degree+1, 0.0, 1.0) + ! Compute bsplines at gauss_points + do k=1,degree+1 + call sll_s_uniform_bsplines_eval_basis(degree,xw_gauss(1,k), bspl(k,:)) + !print*, 'bs', bspl(k,:) + end do + + ! Compute coefs_dofs = int f(x)N_i(x) + do i = 1, self%n_dofs + coef=0.0 + ! loop over support of B spline + do j = 1, degree+1 + ! loop over Gauss points + do k=1, degree+1 + coef = coef + xw_gauss(2,k)*func(self%delta_x*(xw_gauss(1,k) + i + j - 2)) * bspl(k,degree+2-j) + !print*, i,j,k, xw_gauss(2,k), xw_gauss(1,k),f(self%delta_x*(xw_gauss(1,k) + i + j - 2)) + enddo + enddo + ! rescale by cell size + coefs_dofs(i) = coef*self%delta_x + enddo + + end subroutine sll_s_compute_fem_rhs + + !> Compute the L2 projection of a given function f on periodic splines of given degree + subroutine L2projection_1d_fem(self, func, degree, coefs_dofs) + class(sll_t_maxwell_1d_fem) :: self + procedure(sll_i_function_1d_real64) :: funcFFTW.FFTWPlan + sll_int32, intent(in) :: degree + sll_real64, intent(out) :: coefs_dofs(:) ! spline coefficients of projection + ! local variables + sll_int32 :: i + !sll_real64 :: coef + !sll_real64, dimension(2,degree+1) :: xw_gauss + !sll_real64, dimension(degree+1,degree+1) :: bspl + sll_real64, dimension(self%n_dofs) :: eigvals + + ! Compute right-hand-side + call sll_s_compute_fem_rhs(self, func, degree, self%work) + + ! Multiply by inverse mass matrix (! complex numbers stored in real array with fftpack ordering) + eigvals=0.0 + if (degree == self%s_deg_0) then + do i=1,self%n_dofs/2+1 + eigvals(i) = 1.0 / self%eig_mass0(i) + end do + elseif (degree == self%s_deg_0-1) then + do i=1,self%n_dofs/2+1 + eigvals(i) = 1.0 / self%eig_mass1(i) + end do + else + print*, 'degree ', degree, 'not availlable in maxwell_1d_fem object' + endif + + call solve_circulant(self, eigvals, self%work, coefs_dofs) + ! Account for scaling in the mass matrix by dx + coefs_dofs = coefs_dofs/self%delta_x + + end subroutine L2projection_1d_fem + + !> Compute square of the L2norm + function L2norm_squared_1d_fem(self, coefs_dofs, degree) result (r) + class(sll_t_maxwell_1d_fem) :: self !< Maxwell solver object + sll_real64 :: coefs_dofs(:) !< Coefficient for each DoF + sll_int32 :: degree !< Specify the degree of the basis functions + sll_real64 :: r !< Result: squared L2 norm + + ! Multiply coefficients by mass matrix (use diagonalization FFT and mass matrix eigenvalues) + if (degree == self%s_deg_0 ) then + + call solve_circulant(self, self%eig_mass0, coefs_dofs, self%work) + + elseif (degree == self%s_deg_1) then + + call solve_circulant(self, self%eig_mass1, coefs_dofs, self%work) + + end if + ! Multiply by the coefficients from the left (inner product) + r = sum(coefs_dofs*self%work) + ! Scale by delt_x + r = r*self%delta_x + + end function L2norm_squared_1d_fem + + function inner_product_1d_fem( self, coefs1_dofs, coefs2_dofs, degree ) result (r) + class(sll_t_maxwell_1d_fem) :: self !< Maxwell solver object + sll_real64 :: coefs1_dofs(:) !< Coefficient for each DoF + sll_real64 :: coefs2_dofs(:) !< Coefficient for each DoF + sll_int32 :: degree !< Specify the degree of the basis functions + sll_real64 :: r !< Result: squared L2 norm + + ! Multiply coefficients by mass matrix (use diagonalization FFT and mass matrix eigenvalues) + if (degree == self%s_deg_0 ) then + + call solve_circulant(self, self%eig_mass0, coefs2_dofs, self%work) + + elseif (degree == self%s_deg_1) then + + call solve_circulant(self, self%eig_mass1, coefs2_dofs, self%work) + + end if + ! Multiply by the coefficients from the left (inner product) + r = sum(coefs1_dofs*self%work) + ! Scale by delt_x + r = r*self%delta_x + + end function inner_product_1d_fem + + + + subroutine free_1d_fem(self) + class(sll_t_maxwell_1d_fem) :: self + + call sll_s_fft_free( self%plan_fw ) + call sll_s_fft_free( self%plan_bw ) + deallocate(self%mass_0) + deallocate(self%mass_1) + deallocate(self%eig_mass0) + deallocate(self%eig_mass1) + deallocate(self%eig_weak_ampere) + deallocate(self%eig_weak_poisson) + deallocate(self%wsave) + deallocate(self%work) + + end subroutine free_1d_fem + + +end module sll_m_maxwell_1d_fem + +=# diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..cc04dbb --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,177 @@ +using Test +using GEMPIC + +""" + test_maxwell_1d_fem() + + test 1D Maxwell spline finite element solver on a periodic grid + +`L_x` domain dimensions and M is an integer. + +```math +B_z(x,y,t) = \\cos(\\frac{2 M \\pi}{L_x} x) \\cos(\\frac{2 M \\pi}{L_x} t) +``` +```math +E_y(x,y,t) = \\sin(\\frac{2 M \\pi}{L_x} x) \\sin(\\frac{2 M \\pi}{L_x} t) +``` + +""" +function test_maxwell_1d_fem() + + #sll_m_low_level_bsplines, only: sll_s_eval_uniform_periodic_spline_curve + + + eta1_min = .0 + eta1_max = 2π + nc_eta1 = 128 + + Lx = eta1_max - eta1_min + + delta_eta1 = Lx / nc_eta1 + + domain = [eta1_min, eta1_max] + + deg = 3 + + maxwell_1d = Maxwell1DFEM(domain, nc_eta1, deg) + +#= + ! Allocate arrays + SLL_CLEAR_ALLOCATE(ex(1:nc_eta1),error) + SLL_CLEAR_ALLOCATE(ey(1:nc_eta1),error) + SLL_CLEAR_ALLOCATE(bz(1:nc_eta1),error) + + SLL_CLEAR_ALLOCATE(bz_exact(1:nc_eta1),error) + SLL_CLEAR_ALLOCATE(ex_exact(1:nc_eta1),error) + SLL_CLEAR_ALLOCATE(ey_exact(1:nc_eta1),error) + SLL_CLEAR_ALLOCATE(rho(1:nc_eta1),error) + + SLL_CLEAR_ALLOCATE(sval(1:nc_eta1),error) + + + ! Test Poisson + !------------- + ! Set exact solution + do i = 1, nc_eta1 + xi = eta1_min + real(i-1,f64)*delta_eta1 + ex_exact(i) = sin_k(xi)/(2.0_f64*mode*sll_p_pi/Lx) + end do + + call maxwell_1d%compute_rhs_from_function( cos_k, deg, rho) + + call maxwell_1d%compute_e_from_rho( ex, rho ) + + ! Evaluate spline curve at grid points and compute error + ! Ex is a 1-form, i.e. one spline degree lower + call sll_s_eval_uniform_periodic_spline_curve(deg-1, ex, sval) + err_ex = maxval(sval-ex_exact) + print*, 'error Poisson', err_ex + + call sll_s_plot_two_fields_1d('ex',nc_eta1,sval,ex_exact,0,0.0_f64) + + ! Test Ampere + !------------- + ! Set time step + dt = .5_f64 * delta_eta1 + ! Set exact solution + do i = 1, nc_eta1 + xi = eta1_min + real(i-1,f64)*delta_eta1 + ex_exact(i) = -cos_k(xi)*dt + end do + + call maxwell_1d%compute_rhs_from_function(cos_k, deg-1, rho) + ex = 0.0_f64 + call maxwell_1d%compute_E_from_j(dt*rho, 1, ex ) + + ! Evaluate spline curve at grid points and compute error + ! Ex is a 1-form, i.e. one spline degree lower + call sll_s_eval_uniform_periodic_spline_curve(deg-1, ex, sval) + err_ex2 = maxval(sval-ex_exact) + print*, 'error Ampere', err_ex2 + !call sll_plot_two_fields_1d('ex',nc_eta1,sval,ex_exact,0,0.0_f64) + + l2norm = maxwell_1d%l2norm_squared(ex,deg-1) + err_l2norm = l2norm - dt**2*sll_p_pi + print*, 'error l2 norm', err_l2norm + + + ! Test Maxwell on By and Ez + !-------------------------- + ! Set time stepping parameters + time = 0.0_f64 + dt = .5_f64 * delta_eta1 + nstep = 10 + + ! Compute initial fields + ex = 0.0_f64 ! 0-form -> splines of degree deg + call maxwell_1d%L2projection(cos_k, deg-1, bz) ! 0-form -> splines of degree deg-1 + + ! Time loop. Second order Strang splitting + do istep = 1, nstep + call maxwell_1d%compute_b_from_e( 0.5_f64*dt, ey, bz) + call maxwell_1d%compute_e_from_b( dt, bz, ey) + call maxwell_1d%compute_b_from_e( 0.5_f64*dt, ey, bz) + + time = time + dt + + do i = 1, nc_eta1 + xi = eta1_min + real(i-1,f64)*delta_eta1 + ey_exact(i) = sin(mode*2*sll_p_pi*xi/Lx) * sin(mode*2*sll_p_pi*time/Lx) + bz_exact(i) = cos(mode*2*sll_p_pi*xi/Lx) * cos(mode*2*sll_p_pi*time/Lx) + end do + + call sll_s_eval_uniform_periodic_spline_curve(deg, ey, sval) + err_ey = maxval(sval-ey_exact) + call sll_s_eval_uniform_periodic_spline_curve(deg-1, bz, sval) + err_bz = maxval(sval-bz_exact) + + write(*,"(10x,' istep = ',I6)",advance="no") istep + write(*,"(' time = ',g12.3,' sec')",advance="no") time + write(*,"(' erreur L2 = ',2g15.5)") err_ey, err_bz + + call sll_s_plot_two_fields_1d('bz',nc_eta1,sval,bz_exact,istep,time) + + end do ! next time step + + tol = 1.0d-3 + if ((err_bz < tol) .and. (err_ey < tol) .and. (err_ex < tol) .and. (err_ex2 < tol) .and. (err_l2norm < tol)) then + print*,'PASSED' + endif + + call maxwell_1d%free() + + DEALLOCATE(ex) + DEALLOCATE(ey) + DEALLOCATE(bz) + DEALLOCATE(bz_exact) + DEALLOCATE(ex_exact) + DEALLOCATE(ey_exact) + DEALLOCATE(rho) + +contains + + function cos_k(x) + + sll_real64 :: cos_k + sll_real64, intent(in) :: x + + cos_k = cos(mode*2*sll_p_pi*x/Lx) + + end function cos_k + + function sin_k(x) + + sll_real64 :: sin_k + sll_real64, intent(in) :: x + + sin_k = sin(mode*2*sll_p_pi*x/Lx) + + end function sin_k + +=# + + true + +end + +@test test_maxwell_1d_fem()