@@ -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
@@ -151,19 +152,20 @@ pure function output(self) result(state)
151152 ! integrand_tester_object deferred methods
152153 pure function description (self , prefix ) result(desc)
153154 ! < Return informative integrator description.
154- class(integrand_ladvection), intent (in ) :: self ! < Integrator .
155+ class(integrand_ladvection), intent (in ) :: self ! < Integrand .
155156 character (* ), intent (in ), optional :: prefix ! < Prefixing string.
156157 character (len= :), allocatable :: desc ! < Description.
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
163- pure function error (self , t , U0 )
164+ pure function error (self , t , t0 , U0 )
164165 ! < Return error.
165166 class(integrand_ladvection), intent (in ) :: self ! < Integrand.
166167 real (R_P), intent (in ) :: t ! < Time.
168+ real (R_P), intent (in ), optional :: t0 ! < Initial time.
167169 class(integrand_object), intent (in ), optional :: U0 ! < Initial conditions.
168170 real (R_P), allocatable :: error(:) ! < Error.
169171 integer (I_P) :: i ! < Counter.
@@ -174,48 +176,62 @@ pure function error(self, t, U0)
174176 select type (U0)
175177 type is (integrand_ladvection)
176178 do i= 1 , self% Ni
177- ! error = error + (U0%u(i) - self%u(i)) ** 2
178- error = max (error, abs (U0% u(i) - self% u(i)))
179+ error = error + (U0% u(i) - self% u(i)) ** 2
179180 enddo
180181 endselect
181- ! error = sqrt(error)
182+ error = self % Dx * sqrt (error)
182183 endif
183184 endfunction error
184185
185- pure function exact_solution (self , t , U0 ) result(exact)
186+ pure function exact_solution (self , t , t0 , U0 ) result(exact)
186187 ! < Return exact solution.
187188 class(integrand_ladvection), intent (in ) :: self ! < Integrand.
188189 real (R_P), intent (in ) :: t ! < Time.
190+ real (R_P), intent (in ), optional :: t0 ! < Initial time.
189191 class(integrand_object), intent (in ), optional :: U0 ! < Initial conditions.
190192 real (R_P), allocatable :: exact(:) ! < Exact solution.
193+ integer (I_P) :: offset ! < Cells offset.
194+ integer (I_P) :: i ! < Counter.
191195
196+ allocate (exact(1 :self% Ni))
192197 if (present (U0)) then
193198 select type (U0)
194199 type is (integrand_ladvection)
195- 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
196207 endselect
197208 else
198209 exact = self% u(1 :self% Ni) * 0._R_P
199210 endif
200211 endfunction exact_solution
201212
202- 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 )
203214 ! < Export integrand to Tecplot file.
204- class(integrand_ladvection), intent (in ) :: self ! < Advection field.
205- character (* ), intent (in ), optional :: file_name ! < File name.
206- real (R_P), intent (in ), optional :: t ! < Time.
207- character (* ), intent (in ), optional :: scheme ! < Scheme used to integrate integrand.
208- logical , intent (in ), optional :: close_file ! < Flag for closing file.
209- logical , save :: is_open= .false. ! < Flag for checking if file is open.
210- integer (I_P), save :: file_unit ! < File unit.
211- 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.
212227
213228 if (present (close_file)) then
214229 if (close_file .and. is_open) then
215230 close (unit= file_unit)
216231 is_open = .false.
217232 endif
218233 else
234+ with_exact_solution_ = .false. ; if (present (with_exact_solution)) with_exact_solution_ = with_exact_solution
219235 if (present (file_name)) then
220236 if (is_open) close (unit= file_unit)
221237 open (newunit= file_unit, file= trim (adjustl (file_name)))
@@ -227,6 +243,13 @@ subroutine export_tecplot(self, file_name, t, scheme, close_file)
227243 do i= 1 , self% Ni
228244 write (unit= file_unit, fmt= ' (2(' // FR_P// ' ,1X))' ) self% Dx * i - 0.5_R_P * self% Dx, self% u(i)
229245 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
230253 endif
231254 endif
232255 end subroutine export_tecplot
@@ -240,6 +263,8 @@ subroutine initialize(self, Dt)
240263 self% Ni = nint (self% length / self% Dx)
241264
242265 select case (trim (adjustl (self% initial_state)))
266+ case (' sin_wave' )
267+ call self% set_sin_wave_initial_state
243268 case (' square_wave' )
244269 call self% set_square_wave_initial_state
245270 endselect
@@ -252,14 +277,14 @@ subroutine parse_cli(self, cli)
252277
253278 call self% destroy
254279
255- call cli% get(switch= ' --cfl' , val= self% CFL, error= cli% error) ; if (cli% error/= 0 ) stop
256- call cli% get(switch= ' --w-scheme' , val= self% w_scheme, error= cli% error) ; if (cli% error/= 0 ) stop
257- call cli% get(switch= ' --weno-order' , val= self% weno_order, error= cli% error) ; if (cli% error/= 0 ) stop
258- call cli% get(switch= ' --weno-eps' , val= self% weno_eps, error= cli% error) ; if (cli% error/= 0 ) stop
259- call cli% get(switch= ' -a' , val= self% a, error= cli% error) ; if (cli% error/= 0 ) stop
260- call cli% get(switch= ' --length' , val= self% length, error= cli% error) ; if (cli% error/= 0 ) stop
261- call cli% get(switch= ' --Ni' , val= self% Ni, error= cli% error) ; if (cli% error/= 0 ) stop
262- call cli% get(switch= ' -is' , val= self% initial_state, error= cli% error) ; if (cli% error/= 0 ) stop
280+ call cli% get(group = ' linear_advection ' , switch= ' --cfl' , val= self% CFL, error= cli% error) ; if (cli% error/= 0 ) stop
281+ call cli% get(group = ' linear_advection ' , switch= ' --w-scheme' , val= self% w_scheme, error= cli% error) ; if (cli% error/= 0 ) stop
282+ call cli% get(group = ' linear_advection ' , switch= ' --weno-order' , val= self% weno_order, error= cli% error) ; if (cli% error/= 0 ) stop
283+ call cli% get(group = ' linear_advection ' , switch= ' --weno-eps' , val= self% weno_eps, error= cli% error) ; if (cli% error/= 0 ) stop
284+ call cli% get(group = ' linear_advection ' , switch= ' -a' , val= self% a, error= cli% error) ; if (cli% error/= 0 ) stop
285+ call cli% get(group = ' linear_advection ' , switch= ' --length' , val= self% length, error= cli% error) ; if (cli% error/= 0 ) stop
286+ call cli% get(group = ' linear_advection ' , switch= ' --Ni' , val= self% Ni, error= cli% error) ; if (cli% error/= 0 ) stop
287+ call cli% get(group = ' linear_advection ' , switch= ' -is' , val= self% initial_state, error= cli% error) ; if (cli% error/= 0 ) stop
263288
264289 self% Ng = (self% weno_order + 1 ) / 2
265290 self% Dx = self% length / self% Ni
@@ -274,16 +299,19 @@ subroutine set_cli(cli)
274299 ! < Set command line interface.
275300 type (command_line_interface), intent (inout ) :: cli ! < Command line interface handler.
276301
277- call cli% add(switch= ' --w-scheme' , help= ' WENO scheme' , required= .false. , act= ' store' , def= ' reconstructor-JS' , &
278- choices= ' reconstructor-JS,reconstructor-M-JS,reconstructor-M-Z,reconstructor-Z' )
279- call cli% add(switch= ' --weno-order' , help= ' WENO order' , required= .false. , act= ' store' , def= ' 1' )
280- call cli% add(switch= ' --weno-eps' , help= ' WENO epsilon parameter' , required= .false. , act= ' store' , def= ' 0.000001' )
281- call cli% add(switch= ' --cfl' , help= ' CFL value' , required= .false. , act= ' store' , def= ' 0.8' )
282- call cli% add(switch= ' -a' , help= ' advection coefficient' , required= .false. , act= ' store' , def= ' 1.0' )
283- call cli% add(switch= ' --length' , help= ' domain lenth' , required= .false. , act= ' store' , def= ' 1.0' )
284- call cli% add(switch= ' --Ni' , help= ' number finite volumes used' , required= .false. , act= ' store' , def= ' 100' )
285- call cli% add(switch= ' --initial_state' , switch_ab= ' -is' , help= ' initial state' , required= .false. , act= ' store' , &
286- def= ' square_wave' , choices= ' square_wave' )
302+ call cli% add_group(description= ' linear advection test settings' , group= ' linear_advection' )
303+ call cli% add(group= ' linear_advection' , switch= ' --w-scheme' , help= ' WENO scheme' , required= .false. , act= ' store' , &
304+ def= ' reconstructor-JS' , choices= ' reconstructor-JS,reconstructor-M-JS,reconstructor-M-Z,reconstructor-Z' )
305+ call cli% add(group= ' linear_advection' , switch= ' --weno-order' , help= ' WENO order' , required= .false. , act= ' store' , def= ' 1' )
306+ call cli% add(group= ' linear_advection' , switch= ' --weno-eps' , help= ' WENO epsilon parameter' , required= .false. , act= ' store' , &
307+ def= ' 0.000001' )
308+ call cli% add(group= ' linear_advection' , switch= ' --cfl' , help= ' CFL value' , required= .false. , act= ' store' , def= ' 0.8' )
309+ call cli% add(group= ' linear_advection' , switch= ' -a' , help= ' advection coefficient' , required= .false. , act= ' store' , def= ' 1.0' )
310+ call cli% add(group= ' linear_advection' , switch= ' --length' , help= ' domain lenth' , required= .false. , act= ' store' , def= ' 1.0' )
311+ call cli% add(group= ' linear_advection' , switch= ' --Ni' , help= ' number finite volumes used' , required= .false. , act= ' store' , &
312+ def= ' 100' )
313+ call cli% add(group= ' linear_advection' , switch= ' --initial_state' , switch_ab= ' -is' , help= ' initial state' , required= .false. , &
314+ act= ' store' , def= ' sin_wave' , choices= ' sin_wave,square_wave' )
287315 end subroutine set_cli
288316
289317 ! integrand_object deferred methods
@@ -305,7 +333,7 @@ function dU_dt(self, t) result(dState_dt)
305333 do i= 0 , self% Ni
306334 call solve_riemann_problem(state_left= ur(2 , i), state_right= ur(1 , i+1 ), flux= f(i))
307335 enddo
308- allocate (dState_dt(1 - self % Ng :self% Ni+ self % Ng ))
336+ allocate (dState_dt(1 :self% Ni))
309337 do i= 1 , self% Ni
310338 dState_dt(i) = (f(i - 1 ) - f(i)) / self% Dx
311339 enddo
@@ -556,20 +584,34 @@ subroutine reconstruct_interfaces(self, conservative, r_conservative)
556584
557585 subroutine set_square_wave_initial_state (self )
558586 ! < Set initial state as a square wave.
559- class(integrand_ladvection), intent (inout ) :: self ! < Advection field.
560- real (R_P) :: x( 1 :self % ni) ! < Cell center x-abscissa values.
561- 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.
562590
563- 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))
564592 do i= 1 , self% Ni
565- x(i) = self% Dx * i - 0.5_R_P * self% Dx
566- 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
567595 self% u(i) = 0._R_P
568- 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
569597 self% u(i) = 1._R_P
570598 else
571599 self% u(i) = 0._R_P
572600 endif
573601 enddo
574602 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
575617endmodule foodie_test_integrand_ladvection
0 commit comments