@@ -101,6 +101,7 @@ module foodie_test_integrand_ladvection
101101 ! private methods
102102 procedure , pass(self), private :: impose_boundary_conditions ! < Impose boundary conditions.
103103 procedure , pass(self), private :: reconstruct_interfaces ! < Reconstruct interface states.
104+ procedure , pass(self), private :: set_sin_wave_initial_state ! < Set initial state as a sin wave.
104105 procedure , pass(self), private :: set_square_wave_initial_state ! < Set initial state as a square wave.
105106endtype integrand_ladvection
106107
@@ -157,7 +158,7 @@ pure function description(self, prefix) result(desc)
157158 character (len= :), allocatable :: prefix_ ! < Prefixing string, local variable.
158159
159160 prefix_ = ' ' ; if (present (prefix)) prefix_ = prefix
160- desc = prefix// ' linear_advection-Ni_' // trim (strz(self% Ni, 10 ))
161+ desc = prefix// ' linear_advection-' // trim (self % initial_state) // ' - Ni_' // trim (strz(self% Ni, 10 ))
161162 end function description
162163
163164 pure function error (self , t , t0 , U0 )
@@ -175,11 +176,10 @@ pure function error(self, t, t0, U0)
175176 select type (U0)
176177 type is (integrand_ladvection)
177178 do i= 1 , self% Ni
178- ! error = error + (U0%u(i) - self%u(i)) ** 2
179- error = max (error, abs (U0% u(i) - self% u(i)))
179+ error = error + (U0% u(i) - self% u(i)) ** 2
180180 enddo
181181 endselect
182- ! error = sqrt(error)
182+ error = self % Dx * sqrt (error)
183183 endif
184184 endfunction error
185185
@@ -190,34 +190,48 @@ pure function exact_solution(self, t, t0, U0) result(exact)
190190 real (R_P), intent (in ), optional :: t0 ! < Initial time.
191191 class(integrand_object), intent (in ), optional :: U0 ! < Initial conditions.
192192 real (R_P), allocatable :: exact(:) ! < Exact solution.
193+ integer (I_P) :: offset ! < Cells offset.
194+ integer (I_P) :: i ! < Counter.
193195
196+ allocate (exact(1 :self% Ni))
194197 if (present (U0)) then
195198 select type (U0)
196199 type is (integrand_ladvection)
197- exact = U0% u(1 :self% Ni)
200+ offset = nint (mod (self% a * t, self% length) / self% Dx)
201+ do i= 1 , self% Ni - offset
202+ exact(i + offset) = U0% u(i)
203+ enddo
204+ do i= self% Ni - offset + 1 , self% Ni
205+ exact(i - self% Ni + offset) = U0% u(i)
206+ enddo
198207 endselect
199208 else
200209 exact = self% u(1 :self% Ni) * 0._R_P
201210 endif
202211 endfunction exact_solution
203212
204- subroutine export_tecplot (self , file_name , t , scheme , close_file )
213+ subroutine export_tecplot (self , file_name , t , scheme , close_file , with_exact_solution , U0 )
205214 ! < Export integrand to Tecplot file.
206- class(integrand_ladvection), intent (in ) :: self ! < Advection field.
207- character (* ), intent (in ), optional :: file_name ! < File name.
208- real (R_P), intent (in ), optional :: t ! < Time.
209- character (* ), intent (in ), optional :: scheme ! < Scheme used to integrate integrand.
210- logical , intent (in ), optional :: close_file ! < Flag for closing file.
211- logical , save :: is_open= .false. ! < Flag for checking if file is open.
212- integer (I_P), save :: file_unit ! < File unit.
213- integer (I_P) :: i ! < Counter.
215+ class(integrand_ladvection), intent (in ) :: self ! < Advection field.
216+ character (* ), intent (in ), optional :: file_name ! < File name.
217+ real (R_P), intent (in ), optional :: t ! < Time.
218+ character (* ), intent (in ), optional :: scheme ! < Scheme used to integrate integrand.
219+ logical , intent (in ), optional :: close_file ! < Flag for closing file.
220+ logical , intent (in ), optional :: with_exact_solution ! < Flag for export also exact solution.
221+ class(integrand_object), intent (in ), optional :: U0 ! < Initial conditions.
222+ logical :: with_exact_solution_ ! < Flag for export also exact solution, local variable.
223+ logical , save :: is_open= .false. ! < Flag for checking if file is open.
224+ integer (I_P), save :: file_unit ! < File unit.
225+ real (R_P), allocatable :: exact_solution(:) ! < Exact solution.
226+ integer (I_P) :: i ! < Counter.
214227
215228 if (present (close_file)) then
216229 if (close_file .and. is_open) then
217230 close (unit= file_unit)
218231 is_open = .false.
219232 endif
220233 else
234+ with_exact_solution_ = .false. ; if (present (with_exact_solution)) with_exact_solution_ = with_exact_solution
221235 if (present (file_name)) then
222236 if (is_open) close (unit= file_unit)
223237 open (newunit= file_unit, file= trim (adjustl (file_name)))
@@ -229,6 +243,13 @@ subroutine export_tecplot(self, file_name, t, scheme, close_file)
229243 do i= 1 , self% Ni
230244 write (unit= file_unit, fmt= ' (2(' // FR_P// ' ,1X))' ) self% Dx * i - 0.5_R_P * self% Dx, self% u(i)
231245 enddo
246+ if (with_exact_solution_) then
247+ exact_solution = self% exact_solution(t= t, U0= U0)
248+ write (unit= file_unit, fmt= ' (A)' ) ' ZONE T="' // str(t)// ' exact solution"'
249+ do i= 1 , self% Ni
250+ write (unit= file_unit, fmt= ' (2(' // FR_P// ' ,1X))' ) self% Dx * i - 0.5_R_P * self% Dx, exact_solution(i)
251+ enddo
252+ endif
232253 endif
233254 endif
234255 end subroutine export_tecplot
@@ -242,6 +263,8 @@ subroutine initialize(self, Dt)
242263 self% Ni = nint (self% length / self% Dx)
243264
244265 select case (trim (adjustl (self% initial_state)))
266+ case (' sin_wave' )
267+ call self% set_sin_wave_initial_state
245268 case (' square_wave' )
246269 call self% set_square_wave_initial_state
247270 endselect
@@ -288,7 +311,7 @@ subroutine set_cli(cli)
288311 call cli% add(group= ' linear_advection' , switch= ' --Ni' , help= ' number finite volumes used' , required= .false. , act= ' store' , &
289312 def= ' 100' )
290313 call cli% add(group= ' linear_advection' , switch= ' --initial_state' , switch_ab= ' -is' , help= ' initial state' , required= .false. , &
291- act= ' store' , def= ' square_wave ' , choices= ' square_wave' )
314+ act= ' store' , def= ' sin_wave ' , choices= ' sin_wave, square_wave' )
292315 end subroutine set_cli
293316
294317 ! integrand_object deferred methods
@@ -310,7 +333,7 @@ function dU_dt(self, t) result(dState_dt)
310333 do i= 0 , self% Ni
311334 call solve_riemann_problem(state_left= ur(2 , i), state_right= ur(1 , i+1 ), flux= f(i))
312335 enddo
313- allocate (dState_dt(1 - self % Ng :self% Ni+ self % Ng ))
336+ allocate (dState_dt(1 :self% Ni))
314337 do i= 1 , self% Ni
315338 dState_dt(i) = (f(i - 1 ) - f(i)) / self% Dx
316339 enddo
@@ -561,20 +584,34 @@ subroutine reconstruct_interfaces(self, conservative, r_conservative)
561584
562585 subroutine set_square_wave_initial_state (self )
563586 ! < Set initial state as a square wave.
564- class(integrand_ladvection), intent (inout ) :: self ! < Advection field.
565- real (R_P) :: x( 1 :self % ni) ! < Cell center x-abscissa values.
566- integer (I_P) :: i ! < Space counter.
587+ class(integrand_ladvection), intent (inout ) :: self ! < Advection field.
588+ real (R_P) :: x ! < Cell center x-abscissa values.
589+ integer (I_P) :: i ! < Space counter.
567590
568- if (allocated (self% u)) deallocate (self% u) ; allocate (self% u(1 - self % Ng :self% Ni+ self % Ng ))
591+ if (allocated (self% u)) deallocate (self% u) ; allocate (self% u(1 :self% Ni))
569592 do i= 1 , self% Ni
570- x(i) = self% Dx * i - 0.5_R_P * self% Dx
571- if (x(i) < 0.25_R_P ) then
593+ x = self% Dx * i - 0.5_R_P * self% Dx
594+ if (x < 0.25_R_P ) then
572595 self% u(i) = 0._R_P
573- elseif (0.25_R_P <= x(i) .and. x(i) < 0.75_R_P ) then
596+ elseif (0.25_R_P <= x .and. x < 0.75_R_P ) then
574597 self% u(i) = 1._R_P
575598 else
576599 self% u(i) = 0._R_P
577600 endif
578601 enddo
579602 end subroutine set_square_wave_initial_state
603+
604+ subroutine set_sin_wave_initial_state (self )
605+ ! < Set initial state as a sin wave.
606+ class(integrand_ladvection), intent (inout ) :: self ! < Advection field.
607+ real (R_P) :: x ! < Cell center x-abscissa values.
608+ real (R_P), parameter :: pi= 4._R_P * atan (1._R_P ) ! < Pi greek.
609+ integer (I_P) :: i ! < Space counter.
610+
611+ if (allocated (self% u)) deallocate (self% u) ; allocate (self% u(1 :self% Ni))
612+ do i= 1 , self% Ni
613+ x = self% Dx * i - 0.5_R_P * self% Dx
614+ self% u(i) = sin (x * 2 * pi)
615+ enddo
616+ end subroutine set_sin_wave_initial_state
580617endmodule foodie_test_integrand_ladvection
0 commit comments