MFC:Pre_process  v1.0
m_initial_condition.f90
Go to the documentation of this file.
1 !! __ _______________
2 !! / |/ / ____/ ____/
3 !! / /|_/ / /_ / /
4 !! / / / / __/ / /___
5 !! /_/ /_/_/ \____/
6 !!
7 !! This file is part of MFC.
8 !!
9 !! MFC is the legal property of its developers, whose names
10 !! are listed in the copyright file included with this source
11 !! distribution.
12 !!
13 !! MFC is free software: you can redistribute it and/or modify
14 !! it under the terms of the GNU General Public License as published
15 !! by the Free Software Foundation, either version 3 of the license
16 !! or any later version.
17 !!
18 !! MFC is distributed in the hope that it will be useful,
19 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
20 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 !! GNU General Public License for more details.
22 !!
23 !! You should have received a copy of the GNU General Public License
24 !! along with MFC (LICENSE).
25 !! If not, see <http://www.gnu.org/licenses/>.
26 
33 
47 
48 
49  ! Dependencies =============================================================
50  USE m_derived_types ! Definitions of the derived types
51 
52  USE m_global_parameters ! Global parameters for the code
53 
54  USE m_variables_conversion ! Subroutines to change the state variables from
55  ! one form to another
56  ! ==========================================================================
57 
58 
59  IMPLICIT NONE
60 
61 
65  abstract INTERFACE
66 
73  SUBROUTINE s_assign_patch_xxxxx_primitive_variables( patch_id,j,k,l)
74 
75  INTEGER, INTENT(IN) :: patch_id
76  INTEGER, INTENT(IN) :: j,k,l
77 
79 
80  END INTERFACE
81 
82  ! NOTE: The abstract interface allows for the declaration of a pointer to
83  ! a procedure such that the choice of the model equations does not have to
84  ! be queried every time the patch primitive variables are to be assigned in
85  ! a cell in the computational domain.
86  TYPE(scalar_field), ALLOCATABLE, DIMENSION(:) :: q_prim_vf
87 
88  TYPE(scalar_field), ALLOCATABLE, DIMENSION(:) :: q_cons_vf
89 
91 
92  REAL(KIND(0d0)) :: x_centroid, y_centroid, z_centroid
93  REAL(KIND(0d0)) :: length_x, length_y, length_z
94  REAL(KIND(0d0)) :: radius
95  REAL(KIND(0d0)) :: epsilon, beta
96  INTEGER :: smooth_patch_id
97  REAL(KIND(0d0)) :: smooth_coeff
102 
108 
109 
110  REAL(KIND(0d0)) :: a, b, c, d
114 
115 
116  REAL(KIND(0d0)) :: eta
121 
122 
123  INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: patch_id_fp
127 
128  REAL(KIND(0d0)) :: cart_y, cart_z
129  REAL(KIND(0d0)) :: sph_phi
132 
133 
135  POINTER :: s_assign_patch_primitive_variables => null()
139 
140 
141 
142  CONTAINS
143 
158  SUBROUTINE s_assign_patch_mixture_primitive_variables( patch_id,j,k,l)
159 
160 
161  INTEGER, INTENT(IN) :: patch_id
162  INTEGER, INTENT(IN) :: j,k,l
163 
164  REAL(KIND(0d0)) :: rho
165  REAL(KIND(0d0)), DIMENSION(E_idx-mom_idx%beg) :: vel
166  REAL(KIND(0d0)) :: pres
167  REAL(KIND(0d0)) :: gamma
168 
169  INTEGER :: i
170 
171 
172  ! Assigning the mixture primitive variables of a uniform state patch
173  IF(patch_icpp(patch_id)%geometry /= 6) THEN
174 
175  ! Transferring the identity of the smoothing patch
176  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
177 
178  ! Density
179  q_prim_vf(1)%sf(j,k,l) = &
180  eta *patch_icpp( patch_id )%rho &
181  + (1d0 - eta)*patch_icpp(smooth_patch_id)%rho
182 
183  ! Velocity
184  DO i = 1, e_idx-mom_idx%beg
185  q_prim_vf(i+1)%sf(j,k,l) = &
186  1d0 / q_prim_vf(1)%sf(j,k,l) * &
187  ( eta *patch_icpp( patch_id )%rho &
188  *patch_icpp( patch_id )%vel(i) &
189  + (1d0 - eta)*patch_icpp(smooth_patch_id)%rho &
190  *patch_icpp(smooth_patch_id)%vel(i) )
191  END DO
192 
193  ! Specific heat ratio function
194  q_prim_vf(gamma_idx)%sf(j,k,l) = &
195  eta *patch_icpp( patch_id )%gamma &
196  + (1d0 - eta)*patch_icpp(smooth_patch_id)%gamma
197 
198  ! Pressure
199  q_prim_vf(e_idx)%sf(j,k,l) = &
200  1d0 / q_prim_vf(gamma_idx)%sf(j,k,l) * &
201  ( eta *patch_icpp( patch_id )%gamma &
202  *patch_icpp( patch_id )%pres &
203  + (1d0 - eta)*patch_icpp(smooth_patch_id)%gamma &
204  *patch_icpp(smooth_patch_id)%pres )
205 
206  ! Liquid stiffness function
207  q_prim_vf(pi_inf_idx)%sf(j,k,l) = &
208  eta *patch_icpp( patch_id )%pi_inf &
209  + (1d0 - eta)*patch_icpp(smooth_patch_id)%pi_inf
210 
211  ! Assigning mixture primitive variables of isentropic vortex patch
212  ELSE
213 
214  ! Transferring the centroid of the isentropic vortex patch, the
215  ! amplitude of its disturbance and also its domain of influence
216  x_centroid = patch_icpp(patch_id)%x_centroid
217  y_centroid = patch_icpp(patch_id)%y_centroid
218  epsilon = patch_icpp(patch_id)%epsilon
219  beta = patch_icpp(patch_id)%beta
220 
221  ! Reference density, velocity, pressure and specific heat ratio
222  ! function of the isentropic vortex patch
223  rho = patch_icpp(patch_id)%rho
224  vel = patch_icpp(patch_id)%vel
225  pres = patch_icpp(patch_id)%pres
226  gamma = patch_icpp(patch_id)%gamma
227 
228  ! Density
229  q_prim_vf(1)%sf(j,k,0) = &
230  rho*( 1d0 - (rho/pres)*(epsilon/(2d0*pi)) * &
231  (epsilon/(8d0*beta*(gamma+1d0)*pi)) * &
232  exp( 2d0*beta*( 1d0 - (x_cc(j)-x_centroid)**2d0 &
233  - (y_cc(k)-y_centroid)**2d0) ) &
234  )**gamma
235 
236  ! Velocity
237  q_prim_vf(2)%sf(j,k,0) = &
238  vel(1) - (y_cc(k)-y_centroid)*(epsilon/(2d0*pi)) * &
239  exp( beta*( 1d0 - (x_cc(j)-x_centroid)**2d0 &
240  - (y_cc(k)-y_centroid)**2d0) )
241  q_prim_vf(3)%sf(j,k,0) = &
242  vel(2) + (x_cc(j)-x_centroid)*(epsilon/(2d0*pi)) * &
243  exp( beta*( 1d0 - (x_cc(j)-x_centroid)**2d0 &
244  - (y_cc(k)-y_centroid)**2d0) )
245 
246  ! Pressure
247  q_prim_vf(4)%sf(j,k,0) = &
248  pres*( 1d0 - (rho/pres)*(epsilon/(2d0*pi)) * &
249  (epsilon/(8d0*beta*(gamma+1d0)*pi)) * &
250  exp( 2d0*beta*( 1d0 - (x_cc(j)-x_centroid)**2d0 &
251  - (y_cc(k)-y_centroid)**2d0) ) &
252  )**(gamma+1d0)
253 
254  ! Specific heat ratio function
255  q_prim_vf(5)%sf(j,k,0) = gamma
256 
257  ! Liquid stiffness function
258  q_prim_vf(6)%sf(j,k,0) = 0d0
259 
260  END IF
261 
262 
263  ! Updating the patch identities bookkeeping variable
264  IF(1d0 - eta < 1d-16) patch_id_fp(j,k,l) = patch_id
265 
266 
267  END SUBROUTINE s_assign_patch_mixture_primitive_variables ! ------------
268 
276  SUBROUTINE s_assign_patch_species_primitive_variables_bubbles( patch_id,j,k,l)
277 
278  INTEGER, INTENT(IN) :: patch_id
279  INTEGER, INTENT(IN) :: j,k,l
280 
281  ! Density, the specific heat ratio function and the liquid stiffness
282  ! function, respectively, obtained from the combination of primitive
283  ! variables of the current and smoothing patches
284  REAL(KIND(0d0)) :: rho
285  REAL(KIND(0d0)) :: gamma
286  REAL(KIND(0d0)) :: lit_gamma
287  REAL(KIND(0d0)) :: pi_inf
288  REAL(KIND(0d0)) :: orig_rho
289  REAL(KIND(0d0)) :: orig_gamma
290  REAL(KIND(0d0)) :: orig_pi_inf
291 
292  REAL(KIND(0d0)), DIMENSION(sys_size) :: orig_prim_vf
294 
295 
296  INTEGER :: i
297 
298  ! Transferring the identity of the smoothing patch
299  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
300 
301  ! Transferring original primitive variables
302  DO i = 1, sys_size
303  orig_prim_vf(i) = q_prim_vf(i)%sf(j,k,l)
304  END DO
305 
306  IF (mpp_lim .AND. bubbles) THEN
307  !adjust volume fractions, according to modeled gas void fraction
308  alf_sum%sf = 0d0
309  DO i = adv_idx%beg, adv_idx%end - 1
310  alf_sum%sf = alf_sum%sf + q_prim_vf(i)%sf
311  END DO
312 
313  DO i = adv_idx%beg, adv_idx%end-1
314  q_prim_vf(i)%sf = q_prim_vf(i)%sf * (1.d0 - q_prim_vf(alf_idx)%sf) &
315  / alf_sum%sf
316  END DO
317  END IF
318 
319  ! Computing Mixture Variables from Original Primitive Variables
321  q_prim_vf, j,k,l, &
322  orig_rho, &
323  orig_gamma, &
324  orig_pi_inf )
325 
326  ! Computing Mixture Variables of Current Patch =====================
327 
328  ! Volume fraction(s)
329  DO i = adv_idx%beg, adv_idx%end
330  q_prim_vf(i)%sf(j,k,l) = patch_icpp(patch_id)%alpha(i-e_idx)
331  END DO
332 
333  IF (mpp_lim .AND. bubbles) THEN
334  !adjust volume fractions, according to modeled gas void fraction
335  alf_sum%sf = 0d0
336  DO i = adv_idx%beg, adv_idx%end - 1
337  alf_sum%sf = alf_sum%sf + q_prim_vf(i)%sf
338  END DO
339 
340  DO i = adv_idx%beg, adv_idx%end-1
341  q_prim_vf(i)%sf = q_prim_vf(i)%sf * (1.d0 - q_prim_vf(alf_idx)%sf) &
342  / alf_sum%sf
343  END DO
344  END IF
345 
346  ! Partial densities
347  IF (model_eqns .NE. 4) THEN
348  DO i = 1, cont_idx%end
349  q_prim_vf(i)%sf(j,k,l) = patch_icpp(patch_id)%alpha_rho(i)
350  END DO
351  END IF
352 
353  ! Bubbles variables
354  IF (bubbles) THEN
355  DO i = 1,nb
356  q_prim_vf(bub_idx%rs(i))%sf(j,k,l) = r0(i)*patch_icpp(patch_id)%r0
357  q_prim_vf(bub_idx%vs(i))%sf(j,k,l) = v0(i)*patch_icpp(patch_id)%v0
358  END DO
359  END IF
360 
361  ! Density and the specific heat ratio and liquid stiffness functions
363  q_prim_vf, j,k,l, &
364  patch_icpp(patch_id)%rho, &
365  patch_icpp(patch_id)%gamma, &
366  patch_icpp(patch_id)%pi_inf )
367 
368  ! ==================================================================
369 
370  ! Computing Mixture Variables of Smoothing Patch ===================
371 
372  IF (model_eqns .NE. 4) THEN
373  ! Partial densities
374  DO i = 1, cont_idx%end
375  q_prim_vf(i)%sf(j,k,l) = patch_icpp(smooth_patch_id)%alpha_rho(i)
376  END DO
377  END IF
378 
379  ! Volume fraction(s)
380  DO i = adv_idx%beg, adv_idx%end
381  q_prim_vf(i)%sf(j,k,l) = patch_icpp(smooth_patch_id)%alpha(i-e_idx)
382  END DO
383 
384  IF (mpp_lim .AND. bubbles) THEN
385  !adjust volume fractions, according to modeled gas void fraction
386  alf_sum%sf = 0d0
387  DO i = adv_idx%beg, adv_idx%end - 1
388  alf_sum%sf = alf_sum%sf + q_prim_vf(i)%sf
389  END DO
390 
391  DO i = adv_idx%beg, adv_idx%end-1
392  q_prim_vf(i)%sf = q_prim_vf(i)%sf * (1.d0 - q_prim_vf(alf_idx)%sf) &
393  / alf_sum%sf
394  END DO
395  END IF
396 
397  ! Bubbles variables
398  IF (bubbles) THEN
399  DO i = 1,nb
400  q_prim_vf(bub_idx%rs(i))%sf(j,k,l) = r0(i)*patch_icpp(smooth_patch_id)%r0
401  q_prim_vf(bub_idx%vs(i))%sf(j,k,l) = v0(i)*patch_icpp(smooth_patch_id)%v0
402  END DO
403  END IF
404 
405  ! Density and the specific heat ratio and liquid stiffness functions
407  q_prim_vf, j,k,l, &
409  patch_icpp(smooth_patch_id)%gamma, &
410  patch_icpp(smooth_patch_id)%pi_inf )
411 
412  ! ==================================================================
413 
414  ! Pressure
415  q_prim_vf(e_idx)%sf(j,k,l) = &
416  (eta * patch_icpp(patch_id)%pres &
417  + (1d0 - eta) * orig_prim_vf(e_idx))
418 
419  ! Volume fractions \alpha
420  DO i = adv_idx%beg, adv_idx%end
421  q_prim_vf(i)%sf(j,k,l) = &
422  eta *patch_icpp(patch_id)%alpha(i-e_idx) &
423  + (1d0 - eta) *orig_prim_vf(i)
424  END DO
425 
426  IF (mpp_lim .AND. bubbles) THEN
427  !adjust volume fractions, according to modeled gas void fraction
428  alf_sum%sf = 0d0
429  DO i = adv_idx%beg, adv_idx%end - 1
430  alf_sum%sf = alf_sum%sf + q_prim_vf(i)%sf
431  END DO
432 
433  DO i = adv_idx%beg, adv_idx%end-1
434  q_prim_vf(i)%sf = q_prim_vf(i)%sf * (1.d0 - q_prim_vf(alf_idx)%sf) &
435  / alf_sum%sf
436  END DO
437  END IF
438 
439 
440  ! Partial densities \alpha \rho
441  IF (model_eqns .NE. 4) THEN
442  !mixture density is an input
443  DO i = 1, cont_idx%end
444  q_prim_vf(i)%sf(j,k,l) = &
445  eta *patch_icpp(patch_id)%alpha_rho(i) &
446  + (1d0 - eta) * orig_prim_vf(i)
447  END DO
448  ELSE
449  !get mixture density from pressure via Tait EOS
450  pi_inf = fluid_pp(1)%pi_inf
451  gamma = fluid_pp(1)%gamma
452  lit_gamma = (1.d0+gamma)/gamma
453 
454  ! \rho = (( p_l + pi_inf)/( p_ref + pi_inf))**(1/little_gam) * rhoref(1-alf)
455  q_prim_vf(1)%sf(j,k,l) = &
456  (((q_prim_vf(e_idx)%sf(j,k,l) + pi_inf)/(pref + pi_inf))**(1/lit_gamma)) * &
457  rhoref*(1-q_prim_vf(alf_idx)%sf(j,k,l))
458  END IF
459 
460  ! Density and the specific heat ratio and liquid stiffness functions
462  rho,gamma,pi_inf )
463 
464  ! Velocity
465  DO i = 1, e_idx-mom_idx%beg
466  q_prim_vf(i+cont_idx%end)%sf(j,k,l) = &
467  (eta * patch_icpp(patch_id)%vel(i) &
468  + (1d0-eta)*orig_prim_vf(i+cont_idx%end))
469  END DO
470 
471  ! Smoothed bubble variables
472  IF (bubbles) THEN
473  DO i = 1,nb
474  q_prim_vf(bub_idx%rs(i))%sf(j,k,l) = &
475  (eta * r0(i)*patch_icpp(patch_id)%r0 &
476  + (1d0-eta)*orig_prim_vf(bub_idx%rs(i)))
477  q_prim_vf(bub_idx%vs(i))%sf(j,k,l) = &
478  (eta * v0(i)*patch_icpp(patch_id)%v0 &
479  + (1d0-eta)*orig_prim_vf(bub_idx%vs(i)))
480 
481  q_prim_vf(bub_idx%rs(i))%sf(j,k,l) = r0(i)*patch_icpp(patch_id)%r0
482  q_prim_vf(bub_idx%vs(i))%sf(j,k,l) = v0(i)*patch_icpp(patch_id)%v0
483  END DO
484  END IF
485 
486  IF (mpp_lim .AND. bubbles) THEN
487  !adjust volume fractions, according to modeled gas void fraction
488  alf_sum%sf = 0d0
489  DO i = adv_idx%beg, adv_idx%end - 1
490  alf_sum%sf = alf_sum%sf + q_prim_vf(i)%sf
491  END DO
492 
493  DO i = adv_idx%beg, adv_idx%end-1
494  q_prim_vf(i)%sf = q_prim_vf(i)%sf * (1.d0 - q_prim_vf(alf_idx)%sf) &
495  / alf_sum%sf
496  END DO
497  END IF
498 
499  IF (bubbles .AND. (polytropic .NEQV. .true.) ) THEN
500  DO i = 1,nb
501  q_prim_vf(bub_idx%ps(i))%sf(j,k,l) = pb0(i)
502  q_prim_vf(bub_idx%ms(i))%sf(j,k,l) = mass_v0(i)
503  END DO
504  END IF
505 
506 
507  ! Updating the patch identities bookkeeping variable
508  IF(1d0 - eta < 1d-16) patch_id_fp(j,k,l) = patch_id
509 
510  END SUBROUTINE s_assign_patch_species_primitive_variables_bubbles ! ------------
511 
512 
527  SUBROUTINE s_assign_patch_species_primitive_variables( patch_id,j,k,l)
528 
529  INTEGER, INTENT(IN) :: patch_id
530  INTEGER, INTENT(IN) :: j,k,l
531 
532  REAL(KIND(0d0)) :: rho
533  REAL(KIND(0d0)) :: gamma
534  REAL(KIND(0d0)) :: pi_inf
535  REAL(KIND(0d0)) :: orig_rho
536  REAL(KIND(0d0)) :: orig_gamma
537  REAL(KIND(0d0)) :: orig_pi_inf
541 
542  REAL(KIND(0d0)), DIMENSION(sys_size) :: orig_prim_vf
543  ! Vector to hold original values of cell for smoothing purposes
544 
545  INTEGER :: i
546 
547 
548  ! Transferring the identity of the smoothing patch
549  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
550 
551  ! Transferring original primitive variables
552  DO i = 1, sys_size
553  orig_prim_vf(i) = q_prim_vf(i)%sf(j,k,l)
554  END DO
555 
556  ! Computing Mixture Variables from Original Primitive Variables
558  q_prim_vf, j,k,l, &
559  orig_rho, &
560  orig_gamma, &
561  orig_pi_inf )
562 
563  ! Computing Mixture Variables of Current Patch =====================
564 
565  ! Partial densities
566  DO i = 1, cont_idx%end
567  q_prim_vf(i)%sf(j,k,l) = patch_icpp(patch_id)%alpha_rho(i)
568  END DO
569 
570  ! Volume fraction(s)
571  DO i = adv_idx%beg, adv_idx%end
572  q_prim_vf(i)%sf(j,k,l) = patch_icpp(patch_id)%alpha(i-e_idx)
573  END DO
574 
575  ! Density and the specific heat ratio and liquid stiffness functions
577  q_prim_vf, j,k,l, &
578  patch_icpp(patch_id)%rho, &
579  patch_icpp(patch_id)%gamma, &
580  patch_icpp(patch_id)%pi_inf )
581 
582  ! ==================================================================
583 
584  ! Computing Mixture Variables of Smoothing Patch ===================
585 
586  ! Partial densities
587  DO i = 1, cont_idx%end
588  q_prim_vf(i)%sf(j,k,l) = &
589  patch_icpp(smooth_patch_id)%alpha_rho(i)
590  END DO
591 
592  ! Volume fraction(s)
593  DO i = adv_idx%beg, adv_idx%end
594  q_prim_vf(i)%sf(j,k,l) = &
596  END DO
597 
598  ! Density and the specific heat ratio and liquid stiffness functions
600  q_prim_vf, j,k,l, &
602  patch_icpp(smooth_patch_id)%gamma, &
603  patch_icpp(smooth_patch_id)%pi_inf )
604 
605  ! ==================================================================
606 
607  ! Partial densities
608  DO i = 1, cont_idx%end
609  q_prim_vf(i)%sf(j,k,l) = &
610  eta *patch_icpp(patch_id)%alpha_rho(i) &
611  + (1d0 - eta) *orig_prim_vf(i)
612  END DO
613  DO i = adv_idx%beg, adv_idx%end
614  q_prim_vf(i)%sf(j,k,l) = &
615  eta *patch_icpp(patch_id)%alpha(i-e_idx) &
616  + (1d0 - eta) *orig_prim_vf(i)
617  END DO
618 
619  ! Density and the specific heat ratio and liquid stiffness functions
621  rho,gamma,pi_inf )
622 
623  ! Velocity
624  DO i = 1, e_idx-mom_idx%beg
625  q_prim_vf(i+cont_idx%end)%sf(j,k,l) = &
626  (eta * patch_icpp(patch_id)%vel(i) &
627  + (1d0-eta)*orig_prim_vf(i+cont_idx%end))
628  END DO
629 
630  ! Pressure
631  q_prim_vf(e_idx)%sf(j,k,l) = &
632  (eta * patch_icpp(patch_id)%pres &
633  + (1d0 - eta) * orig_prim_vf(e_idx))
634 
635  ! Set partial pressures to mixture pressure
636  IF(model_eqns == 3) THEN
638  q_prim_vf(i)%sf(j,k,l) = q_prim_vf(e_idx)%sf(j,k,l)
639  END DO
640  ENDIF
641 
642  ! Updating the patch identities bookkeeping variable
643  IF(1d0 - eta < 1d-16) patch_id_fp(j,k,l) = patch_id
644 
645 
646  END SUBROUTINE s_assign_patch_species_primitive_variables ! ------------
647 
650  SUBROUTINE s_initialize_initial_condition_module() ! -------------------
651 
652  INTEGER :: i
653 
654 
655  ! Allocating the primitive and conservative variables
656  ALLOCATE(q_prim_vf(1:sys_size))
657  ALLOCATE(q_cons_vf(1:sys_size))
658 
659  DO i = 1, sys_size
660  ALLOCATE(q_prim_vf(i)%sf(0:m,0:n,0:p))
661  ALLOCATE(q_cons_vf(i)%sf(0:m,0:n,0:p))
662  END DO
663  ALLOCATE( alf_sum%sf(0:m,0:n,0:p) )
664 
665 
666  ! Allocating the patch identities bookkeeping variable
667  ALLOCATE(patch_id_fp(0:m,0:n,0:p))
668 
669 
670  ! Setting default values for conservative and primitive variables so
671  ! that in the case that the initial condition is wrongly laid out on
672  ! the grid the simulation component will catch the problem on start-
673  ! up. The conservative variables do not need to be similarly treated
674  ! since they are computed directly from the primitive variables.
675  DO i = 1, sys_size
676  q_cons_vf(i)%sf = dflt_real
677  q_prim_vf(i)%sf = dflt_real
678  END DO
679 
680 
681  ! Setting default values for patch identities bookkeeping variable.
682  ! This is necessary to avoid any confusion in the assessment of the
683  ! extent of application that the overwrite permissions give a patch
684  ! when it is being applied in the domain.
685  patch_id_fp = 0
686 
687 
688  ! Depending on multicomponent flow model, the appropriate procedure
689  ! for assignment of the patch mixture or species primitive variables
690  ! to a cell in the domain is targeted by the procedure pointer
691 
692  IF(model_eqns == 1) THEN ! Gamma/pi_inf model
695  ELSE IF ( bubbles ) THEN
698  ELSE ! Volume fraction model
701  END IF
702 
703 
704  END SUBROUTINE s_initialize_initial_condition_module ! -----------------
705 
706 
707 
714  SUBROUTINE s_generate_initial_condition() ! ----------------------------
715 
716 
717  INTEGER :: i
718 
719 
720  ! Converting the conservative variables to the primitive ones given
721  ! preexisting initial condition data files were read in on start-up
722  IF(old_ic) THEN
724  q_prim_vf )
725  END IF
726 
727 
728  ! 3D Patch Geometries =============================================
729  IF(p > 0) THEN
730 
731  DO i = 1, num_patches
732 
733  ! Spherical patch
734  IF(patch_icpp(i)%geometry == 8) THEN
735  CALL s_sphere(i)
736 
737  ! Cuboidal patch
738  ELSEIF(patch_icpp(i)%geometry == 9) THEN
739  CALL s_cuboid(i)
740 
741  ! Cylindrical patch
742  ELSEIF(patch_icpp(i)%geometry == 10) THEN
743  CALL s_cylinder(i)
744 
745  ! Swept plane patch
746  ELSEIF(patch_icpp(i)%geometry == 11) THEN
747  CALL s_sweep_plane(i)
748 
749  ! Ellipsoidal patch
750  ELSEIF(patch_icpp(i)%geometry == 12) THEN
751  CALL s_ellipsoid(i)
752 
753  ! Analytical function patch for testing purposes
754  ELSEIF(patch_icpp(i)%geometry == 13) THEN
755  CALL s_3d_analytical(i)
756 
757  ! Spherical harmonic patch
758  ELSEIF(patch_icpp(i)%geometry == 14) THEN
759  CALL s_spherical_harmonic(i)
760 
761  ! 3D Modified circular patch
762  ELSEIF(patch_icpp(i)%geometry == 19) THEN
763  CALL s_3dvarcircle(i)
764 
765  ENDIF
766 
767  END DO
768 
769  ! ==================================================================
770 
771 
772  ! 2D Patch Geometries ==============================================
773  ELSEIF(n > 0) THEN
774 
775  DO i = 1, num_patches
776 
777  ! Circular patch
778  IF(patch_icpp(i)%geometry == 2) THEN
779  CALL s_circle(i)
780 
781  ! Rectangular patch
782  ELSEIF(patch_icpp(i)%geometry == 3) THEN
783  CALL s_rectangle(i)
784 
785  ! Swept line patch
786  ELSEIF(patch_icpp(i)%geometry == 4) THEN
787  CALL s_sweep_line(i)
788 
789  ! Elliptical patch
790  ELSEIF(patch_icpp(i)%geometry == 5) THEN
791  CALL s_ellipse(i)
792 
793  ! Isentropic vortex patch
794  ELSEIF(patch_icpp(i)%geometry == 6) THEN
795  CALL s_isentropic_vortex(i)
796 
797  ! Analytical function patch for testing purposes
798  ELSEIF(patch_icpp(i)%geometry == 7) THEN
799  CALL s_2d_analytical(i)
800 
801  ! Spiral patch
802  ELSEIF(patch_icpp(i)%geometry == 17) THEN
803  CALL s_spiral(i)
804 
805  ! Modified circular patch
806  ELSEIF(patch_icpp(i)%geometry == 18) THEN
807  CALL s_varcircle(i)
808 
809 
810  END IF
811 
812  END DO
813 
814  ! ==================================================================
815 
816 
817  ! 1D Patch Geometries ==============================================
818  ELSE
819 
820  DO i = 1, num_patches
821 
822  ! Line segment patch
823  IF (patch_icpp(i)%geometry == 1) THEN
824  CALL s_line_segment(i)
825 
826  ! 1d analytical
827  ELSEIF(patch_icpp(i)%geometry == 15) THEN
828  CALL s_1d_analytical(i)
829 
830  ! 1d bubble screen with sinusoidal pressure pulse
831  ELSEIF(patch_icpp(i)%geometry == 16) THEN
832  CALL s_1d_bubble_pulse(i)
833  END IF
834 
835  END DO
836 
837  END IF
838  ! ==================================================================
839 
841  IF (perturb_sph) CALL s_perturb_sphere()
842 
843  ! Converting the primitive variables to the conservative ones
845  q_cons_vf )
846 
847 
848  END SUBROUTINE s_generate_initial_condition ! --------------------------
849 
850 
851  SUBROUTINE s_convert_cylindrical_to_cartesian_coord(cyl_y,cyl_z)
853  REAL(KIND(0d0)), INTENT(IN) :: cyl_y, cyl_z
854 
855  cart_y = cyl_y*sin(cyl_z)
856  cart_z = cyl_y*cos(cyl_z)
857 
858  END SUBROUTINE s_convert_cylindrical_to_cartesian_coord ! --------------
859 
860  SUBROUTINE s_convert_cylindrical_to_spherical_coord(cyl_x,cyl_y)
862  REAL(KIND(0d0)), INTENT(IN) :: cyl_x, cyl_y
863 
864  sph_phi = atan(cyl_y/cyl_x)
865 
866  END SUBROUTINE s_convert_cylindrical_to_spherical_coord ! --------------
867 
868 
869  SUBROUTINE s_perturb_sphere() ! ----------------------------------------
871  INTEGER :: i,j,k,l
872 
873  REAL(KIND(0d0)) :: perturb_alpha
874  REAL(KIND(0d0)) :: alpha_unadv
875  REAL(KIND(0d0)) :: rand_real
876  CALL random_seed()
877 
878  DO k = 0, p
879  DO j = 0, n
880  DO i = 0, m
881  CALL random_number(rand_real)
882  IF ((perturb_sph_fluid == num_fluids) .AND. (adv_alphan .NEQV. .true.)) THEN
883  perturb_alpha = 1d0
884  DO l = adv_idx%beg, adv_idx%end
885  perturb_alpha = perturb_alpha - q_prim_vf(l)%sf(i,j,k)
886  END DO
887  ELSE
888  perturb_alpha = q_prim_vf(e_idx+perturb_sph_fluid)%sf(i,j,k)
889  END IF
890 
891  ! Perturb partial density fields to match perturbed volume fraction fields
892 ! IF ((perturb_alpha >= 25d-2) .AND. (perturb_alpha <= 75d-2)) THEN
893  IF ((perturb_alpha /= 0d0) .AND. (perturb_alpha /= 1d0)) THEN
894  IF (adv_alphan .NEQV. .true.) THEN
895  ! Find new unadvected volume fraction
896  alpha_unadv = 1d0
897  DO l = adv_idx%beg, adv_idx%end
898  alpha_unadv = alpha_unadv - q_prim_vf(l)%sf(i,j,k)
899  END DO
900  ! Derive new partial densities
901  DO l = 1, num_fluids-1
902  q_prim_vf(l)%sf(i,j,k) = q_prim_vf(e_idx + l)%sf(i,j,k)*fluid_rho(l)
903  END DO
904  q_prim_vf(num_fluids)%sf(i,j,k) = alpha_unadv*fluid_rho(num_fluids)
905  ELSE
906  ! Derive new partial densities
907  DO l = 1, num_fluids
908  q_prim_vf(l)%sf(i,j,k) = q_prim_vf(e_idx + l)%sf(i,j,k)*fluid_rho(l)
909  END DO
910  END IF
911  END IF
912  END DO
913  END DO
914  END DO
915 
916  END SUBROUTINE s_perturb_sphere ! --------------------------------------
917 
918 
919 
920 
921  SUBROUTINE s_perturb_surrounding_flow() ! ------------------------------
923  INTEGER :: i,j,k,l
924 
925  REAL(KIND(0d0)) :: perturb_alpha
926  REAL(KIND(0d0)) :: rand_real
927  CALL random_seed()
928 
929  ! Perturb partial density or velocity of surrounding flow by some random small amount of noise
930  DO k = 0, p
931  DO j = 0, n
932  DO i = 0, m
933  IF ((perturb_flow_fluid == num_fluids) .AND. (adv_alphan .NEQV. .true.)) THEN
934  perturb_alpha = 1d0
935  DO l = adv_idx%beg, adv_idx%end
936  perturb_alpha = perturb_alpha - q_prim_vf(l)%sf(i,j,k)
937  END DO
938  ELSE
939  perturb_alpha = q_prim_vf(e_idx+perturb_flow_fluid)%sf(i,j,k)
940  END IF
941  IF (perturb_alpha == 1d0) THEN
942  ! Perturb partial density
943 ! CALL RANDOM_NUMBER(rand_real)
944 ! rand_real = rand_real / 1d2 / 1d3
945 ! q_prim_vf(perturb_flow_fluid)%sf(i,j,k) = q_prim_vf(perturb_flow_fluid)%sf(i,j,k) + rand_real
946  ! Perturb velocity
947  DO l = mom_idx%beg+1, mom_idx%end
948  CALL random_number(rand_real)
949  rand_real = rand_real / 1d1 / 145d1
950  q_prim_vf(l)%sf(i,j,k) = q_prim_vf(l)%sf(i,j,k) + rand_real
951  END DO
952  END IF
953  END DO
954  END DO
955  END DO
956 
957  END SUBROUTINE s_perturb_surrounding_flow ! ----------------------------
958 
959 
967  SUBROUTINE s_line_segment(patch_id) ! ----------------------------------
969  INTEGER, INTENT(IN) :: patch_id
970 
971  REAL(KIND(0d0)) :: pi_inf, gamma, lit_gamma
972 
973  INTEGER :: i,j
974 
975  pi_inf = fluid_pp(1)%pi_inf
976  gamma = fluid_pp(1)%gamma
977  lit_gamma = (1d0+gamma)/gamma
978 
979  ! Transferring the line segment's centroid and length information
980  x_centroid = patch_icpp(patch_id)%x_centroid
981  length_x = patch_icpp(patch_id)%length_x
982 
983 
984  ! Computing the beginning and end x-coordinates of the line segment
985  ! based on its centroid and length
986  x_boundary%beg = x_centroid - 0.5d0*length_x
987  x_boundary%end = x_centroid + 0.5d0*length_x
988 
989 
990  ! Since the line segment patch does not allow for its boundaries to
991  ! be smoothed out, the pseudo volume fraction is set to 1 to ensure
992  ! that only the current patch contributes to the fluid state in the
993  ! cells that this patch covers.
994  eta = 1d0
995 
996 
997  ! Checking whether the line segment covers a particular cell in the
998  ! domain and verifying whether the current patch has the permission
999  ! to write to that cell. If both queries check out, the primitive
1000  ! variables of the current patch are assigned to this cell.
1001  DO i = 0, m
1002  IF( x_boundary%beg <= x_cc(i) .AND. &
1003  x_boundary%end >= x_cc(i) .AND. &
1004  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,0,0)) ) THEN
1005 
1006  CALL s_assign_patch_primitive_variables(patch_id, i,0,0)
1007 
1008  !IF ( (q_prim_vf(1)%sf(i,0,0) < 1.e-12) .AND. (model_eqns .NE. 4)) THEN
1009  ! !zero density, reassign according to Tait EOS
1010  ! q_prim_vf(1)%sf(i,0,0) = &
1011  ! (((q_prim_vf(E_idx)%sf(i,0,0) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma)) * &
1012  ! rhoref*(1d0-q_prim_vf(alf_idx)%sf(i,0,0))
1013  !END IF
1014  END IF
1015  END DO
1016 
1017  END SUBROUTINE s_line_segment ! ----------------------------------------
1018 
1024  SUBROUTINE s_spiral(patch_id) ! ----------------------------------------
1026  INTEGER, INTENT(IN) :: patch_id
1027 
1028  INTEGER :: i,j,k
1029  REAL(KIND(0d0)) :: th, thickness, nturns, mya
1030  REAL(KIND(0d0)) :: spiral_x_min, spiral_x_max, spiral_y_min, spiral_y_max
1031 
1032  ! Transferring the circular patch's radius, centroid, smearing patch
1033  ! identity and smearing coefficient information
1034  x_centroid = patch_icpp(patch_id)%x_centroid
1035  y_centroid = patch_icpp(patch_id)%y_centroid
1036  mya = patch_icpp(patch_id)%radius
1037  thickness = patch_icpp(patch_id)%length_x
1038  nturns = patch_icpp(patch_id)%length_y
1039 
1040 !
1041  logic_grid = 0
1042  DO k = 0,int(m*91*nturns)
1043  th = k/REAL(int(m*91d0*nturns)) * nturns*2.d0*pi
1044 
1045  spiral_x_min = minval( (/ f_r(th,0.0d0, mya)*cos(th), &
1046  f_r(th,thickness,mya)*cos(th) /) )
1047  spiral_y_min = minval( (/ f_r(th,0.0d0, mya)*sin(th), &
1048  f_r(th,thickness,mya)*sin(th) /) )
1049 
1050  spiral_x_max = maxval( (/ f_r(th,0.0d0, mya)*cos(th), &
1051  f_r(th,thickness,mya)*cos(th) /) )
1052  spiral_y_max = maxval( (/ f_r(th,0.0d0, mya)*sin(th), &
1053  f_r(th,thickness,mya)*sin(th) /) )
1054 
1055  DO j = 0,n; DO i = 0,m;
1056  if( (x_cc(i) > spiral_x_min) .AND. (x_cc(i) < spiral_x_max) .AND. &
1057  (y_cc(j) > spiral_y_min) .AND. (y_cc(j) < spiral_y_max)) THEN
1058  logic_grid(i,j,0) = 1
1059  END IF
1060  END DO; END do
1061  END DO
1062 
1063  DO j = 0, n
1064  DO i = 0, m
1065  IF ( (logic_grid(i,j,0) == 1) ) THEN
1066  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1067  END IF
1068  END DO
1069  END DO
1070 
1071  END SUBROUTINE s_spiral ! ----------------------------------------------
1072 
1073 
1078  FUNCTION f_r(myth,offset,a)
1079  REAL(KIND(0d0)), INTENT(IN) :: myth, offset,a
1080  REAL(KIND(0d0)) :: b
1081  REAL(KIND(0d0)) :: f_r
1082 
1083  !r(th) = a + b*th
1084 
1085  b = 2.d0*a/(2.d0*pi)
1086  f_r = a + b*myth + offset
1087  END FUNCTION f_r
1088 
1095  SUBROUTINE s_circle(patch_id) ! ----------------------------------------
1097  INTEGER, INTENT(IN) :: patch_id
1098 
1099  INTEGER :: i,j
1100 
1101  ! Transferring the circular patch's radius, centroid, smearing patch
1102  ! identity and smearing coefficient information
1103  x_centroid = patch_icpp(patch_id)%x_centroid
1104  y_centroid = patch_icpp(patch_id)%y_centroid
1105  radius = patch_icpp(patch_id)%radius
1106  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
1107  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
1108 
1109 
1110  ! Initializing the pseudo volume fraction value to 1. The value will
1111  ! be modified as the patch is laid out on the grid, but only in the
1112  ! case that smoothing of the circular patch's boundary is enabled.
1113  eta = 1d0
1114 
1115 
1116  ! Checking whether the circle covers a particular cell in the domain
1117  ! and verifying whether the current patch has permission to write to
1118  ! that cell. If both queries check out, the primitive variables of
1119  ! the current patch are assigned to this cell.
1120  DO j = 0, n
1121  DO i = 0, m
1122 
1123  IF(patch_icpp(patch_id)%smoothen) THEN
1124 
1125  eta = tanh( smooth_coeff / min(dx,dy) * &
1126  ( sqrt( (x_cc(i) - x_centroid)**2d0 &
1127  + (y_cc(j) - y_centroid)**2d0 ) &
1128  - radius ) ) * (-0.5d0) + 0.5d0
1129 
1130  END IF
1131 
1132  IF(( (x_cc(i) - x_centroid)**2d0 &
1133  + (y_cc(j) - y_centroid)**2d0 <= radius**2d0 &
1134  .AND. &
1135  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0))) &
1136  .OR. &
1137  patch_id_fp(i,j,0) == smooth_patch_id ) &
1138  THEN
1139 
1140  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1141 
1142  END IF
1143 
1144  END DO
1145  END DO
1146 
1147 
1148  END SUBROUTINE s_circle ! ----------------------------------------------
1149 
1150 
1154  SUBROUTINE s_varcircle(patch_id) ! ----------------------------------------
1156 
1157  ! Patch identifier
1158  INTEGER, INTENT(IN) :: patch_id
1159 
1160  ! Generic loop iterators
1161  INTEGER :: i,j
1162 
1163  REAL(KIND(0d0)) :: myr, thickness
1164 
1165  ! Transferring the circular patch's radius, centroid, smearing patch
1166  ! identity and smearing coefficient information
1167  x_centroid = patch_icpp(patch_id)%x_centroid
1168  y_centroid = patch_icpp(patch_id)%y_centroid
1169  radius = patch_icpp(patch_id)%radius
1170  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
1171  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
1172  thickness = patch_icpp(patch_id)%epsilon
1173 
1174 
1175  ! Initializing the pseudo volume fraction value to 1. The value will
1176  ! be modified as the patch is laid out on the grid, but only in the
1177  ! case that smoothing of the circular patch's boundary is enabled.
1178  eta = 1d0
1179 
1180 
1181  ! Checking whether the circle covers a particular cell in the domain
1182  ! and verifying whether the current patch has permission to write to
1183  ! that cell. If both queries check out, the primitive variables of
1184  ! the current patch are assigned to this cell.
1185  DO j = 0, n
1186  DO i = 0, m
1187  myr = dsqrt( (x_cc(i) - x_centroid)**2d0 &
1188  + (y_cc(j) - y_centroid)**2d0 )
1189 
1190  IF( myr <= radius+thickness/2.d0 .AND. &
1191  myr >= radius-thickness/2.d0 .AND. &
1192  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0)) ) THEN
1193 
1194  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1195 
1196  q_prim_vf(alf_idx)%sf(i,j,0) = patch_icpp(patch_id)%alpha(1) * &
1197  dexp( -0.5d0 * ((myr - radius)**2.d0) / (thickness/3.d0)**2.d0 )
1198  END IF
1199 
1200  END DO
1201  END DO
1202 
1203 
1204  END SUBROUTINE s_varcircle ! ----------------------------------------------
1205 
1206 
1207  SUBROUTINE s_3dvarcircle(patch_id) ! ----------------------------------------
1209  ! Patch identifier
1210  INTEGER, INTENT(IN) :: patch_id
1211 
1212  ! Generic loop iterators
1213  INTEGER :: i,j,k
1214 
1215  REAL(KIND(0d0)) :: myr, thickness
1216 
1217  ! Transferring the circular patch's radius, centroid, smearing patch
1218  ! identity and smearing coefficient information
1219  x_centroid = patch_icpp(patch_id)%x_centroid
1220  y_centroid = patch_icpp(patch_id)%y_centroid
1221  z_centroid = patch_icpp(patch_id)%z_centroid
1222  length_z = patch_icpp(patch_id)%length_z
1223  radius = patch_icpp(patch_id)%radius
1224  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
1225  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
1226  thickness = patch_icpp(patch_id)%epsilon
1227 
1228 
1229  ! Initializing the pseudo volume fraction value to 1. The value will
1230  ! be modified as the patch is laid out on the grid, but only in the
1231  ! case that smoothing of the circular patch's boundary is enabled.
1232  eta = 1d0
1233 
1234  ! write for all z
1235 
1236  ! Checking whether the circle covers a particular cell in the domain
1237  ! and verifying whether the current patch has permission to write to
1238  ! that cell. If both queries check out, the primitive variables of
1239  ! the current patch are assigned to this cell.
1240  DO k =0,p
1241  DO j = 0, n
1242  DO i = 0, m
1243  myr = dsqrt( (x_cc(i) - x_centroid)**2d0 &
1244  + (y_cc(j) - y_centroid)**2d0 )
1245 
1246  IF( myr <= radius+thickness/2.d0 .AND. &
1247  myr >= radius-thickness/2.d0 .AND. &
1248  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k)) ) THEN
1249 
1250  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
1251 
1252  q_prim_vf(alf_idx)%sf(i,j,k) = patch_icpp(patch_id)%alpha(1) * &
1253  dexp( -0.5d0 * ((myr - radius)**2.d0) / (thickness/3.d0)**2.d0 )
1254  END IF
1255 
1256  END DO
1257  END DO
1258  END DO
1259 
1260 
1261  END SUBROUTINE s_3dvarcircle ! ----------------------------------------------
1262 
1268  SUBROUTINE s_ellipse(patch_id) ! ---------------------------------------
1270  INTEGER, INTENT(IN) :: patch_id
1271 
1272  INTEGER :: i,j
1273 
1274  ! Transferring the elliptical patch's radii, centroid, smearing
1275  ! patch identity, and smearing coefficient information
1276  x_centroid = patch_icpp(patch_id)%x_centroid
1277  y_centroid = patch_icpp(patch_id)%y_centroid
1278  a = patch_icpp(patch_id)%radii(1)
1279  b = patch_icpp(patch_id)%radii(2)
1280  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
1281  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
1282 
1283  ! Initializing the pseudo volume fraction value to 1. The value
1284  ! be modified as the patch is laid out on the grid, but only in
1285  ! the case that smoothing of the elliptical patch's boundary is
1286  ! enabled.
1287  eta = 1d0
1288 
1289  ! Checking whether the ellipse covers a particular cell in the
1290  ! domain and verifying whether the current patch has permission
1291  ! to write to that cell. If both queries check out, the primitive
1292  ! variables of the current patch are assigned to this cell.
1293  DO j = 0, n
1294  DO i = 0, m
1295 
1296  IF(patch_icpp(patch_id)%smoothen) THEN
1297  eta = tanh(smooth_coeff / min(dx,dy) * &
1298  ( sqrt( ((x_cc(i) - x_centroid) / a)**2d0 + &
1299  ((y_cc(j) - y_centroid) / b)**2d0 ) &
1300  - 1d0 )) * (-0.5d0) + 0.5d0
1301  END IF
1302 
1303  IF (( ((x_cc(i) - x_centroid) / a)**2d0 + &
1304  ((y_cc(j) - y_centroid) / b)**2d0 <= 1d0 &
1305  .AND. &
1306  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0))) &
1307  .OR. &
1308  patch_id_fp(i,j,0) == smooth_patch_id ) &
1309  THEN
1310 
1311  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1312  END IF
1313  END DO
1314  END DO
1315 
1316  END SUBROUTINE s_ellipse ! ---------------------------------------------
1317 
1318 
1319 
1325  SUBROUTINE s_ellipsoid(patch_id) ! -------------------------------------
1327 
1328  ! Patch identifier
1329  INTEGER, INTENT(IN) :: patch_id
1330 
1331  ! Generic loop iterators
1332  INTEGER :: i,j,k
1333 
1334  ! Transferring the ellipsoidal patch's radii, centroid, smearing
1335  ! patch identity, and smearing coefficient information
1336  x_centroid = patch_icpp(patch_id)%x_centroid
1337  y_centroid = patch_icpp(patch_id)%y_centroid
1338  z_centroid = patch_icpp(patch_id)%z_centroid
1339  a = patch_icpp(patch_id)%radii(1)
1340  b = patch_icpp(patch_id)%radii(2)
1341  c = patch_icpp(patch_id)%radii(3)
1342  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
1343  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
1344 
1345  ! Initializing the pseudo volume fraction value to 1. The value
1346  ! be modified as the patch is laid out on the grid, but only in
1347  ! the case that smoothing of the ellipsoidal patch's boundary is
1348  ! enabled.
1349  eta = 1d0
1350 
1351  ! Checking whether the ellipsoid covers a particular cell in the
1352  ! domain and verifying whether the current patch has permission
1353  ! to write to that cell. If both queries check out, the primitive
1354  ! variables of the current patch are assigned to this cell.
1355  DO k = 0, p
1356  DO j = 0, n
1357  DO i = 0, m
1358 
1359  IF (grid_geometry == 3) THEN
1360  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
1361  ELSE
1362  cart_y = y_cc(j)
1363  cart_z = z_cc(k)
1364  END IF
1365 
1366  IF(patch_icpp(patch_id)%smoothen) THEN
1367  eta = tanh(smooth_coeff / min(dx,dy,dz) * &
1368  ( sqrt( ((x_cc(i) - x_centroid) / a)**2d0 + &
1369  (( cart_y - y_centroid) / b)**2d0 + &
1370  (( cart_z - z_centroid) / c)**2d0 ) &
1371  - 1d0 )) * (-0.5d0) + 0.5d0
1372  END IF
1373 
1374  IF (( (( x_cc(i) - x_centroid) / a)**2d0 + &
1375  (( cart_y - y_centroid) / b)**2d0 + &
1376  (( cart_z - z_centroid) / c)**2d0 <= 1d0 &
1377  .AND. &
1378  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k))) &
1379  .OR. &
1380  patch_id_fp(i,j,k) == smooth_patch_id ) &
1381  THEN
1382 
1383  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
1384  END IF
1385  END DO
1386  END DO
1387  END DO
1388 
1389  END SUBROUTINE s_ellipsoid ! -------------------------------------------
1390 
1391 
1401  SUBROUTINE s_rectangle(patch_id) ! -------------------------------------
1403  INTEGER, INTENT(IN) :: patch_id
1404 
1405  REAL(KIND(0d0)) :: pi_inf, gamma, lit_gamma
1406 
1407  INTEGER :: i,j
1408 
1409  pi_inf = fluid_pp(1)%pi_inf
1410  gamma = fluid_pp(1)%gamma
1411  lit_gamma = (1d0+gamma)/gamma
1412 
1413 
1414 
1415  ! Transferring the rectangle's centroid and length information
1416  x_centroid = patch_icpp(patch_id)%x_centroid
1417  y_centroid = patch_icpp(patch_id)%y_centroid
1418  length_x = patch_icpp(patch_id)%length_x
1419  length_y = patch_icpp(patch_id)%length_y
1420 
1421 
1422  ! Computing the beginning and the end x- and y-coordinates of the
1423  ! rectangle based on its centroid and lengths
1424  x_boundary%beg = x_centroid - 0.5d0*length_x
1425  x_boundary%end = x_centroid + 0.5d0*length_x
1426  y_boundary%beg = y_centroid - 0.5d0*length_y
1427  y_boundary%end = y_centroid + 0.5d0*length_y
1428 
1429 
1430  ! Since the rectangular patch does not allow for its boundaries to
1431  ! be smoothed out, the pseudo volume fraction is set to 1 to ensure
1432  ! that only the current patch contributes to the fluid state in the
1433  ! cells that this patch covers.
1434  eta = 1d0
1435 
1436 
1437  ! Checking whether the rectangle covers a particular cell in the
1438  ! domain and verifying whether the current patch has the permission
1439  ! to write to that cell. If both queries check out, the primitive
1440  ! variables of the current patch are assigned to this cell.
1441  DO j = 0, n
1442  DO i = 0, m
1443  IF( x_boundary%beg <= x_cc(i) .AND. &
1444  x_boundary%end >= x_cc(i) .AND. &
1445  y_boundary%beg <= y_cc(j) .AND. &
1446  y_boundary%end >= y_cc(j) &
1447  .AND. &
1448  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0)) ) &
1449  THEN
1450 
1451  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1452 
1453  IF ( (q_prim_vf(1)%sf(i,j,0) < 1.e-10) .AND. (model_eqns .NE. 4)) THEN
1454  !zero density, reassign according to Tait EOS
1455  q_prim_vf(1)%sf(i,j,0) = &
1456  (((q_prim_vf(e_idx)%sf(i,j,0) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma)) * &
1457  rhoref*(1d0-q_prim_vf(alf_idx)%sf(i,j,0))
1458  END IF
1459  END IF
1460  END DO
1461  END DO
1462 
1463 
1464  END SUBROUTINE s_rectangle ! -------------------------------------------
1465 
1466 
1475  SUBROUTINE s_sweep_line(patch_id) ! ------------------------------------
1477  INTEGER, INTENT(IN) :: patch_id
1478 
1479  INTEGER :: i,j
1480 
1481 
1482  ! Transferring the centroid information of the line to be swept
1483  x_centroid = patch_icpp(patch_id)%x_centroid
1484  y_centroid = patch_icpp(patch_id)%y_centroid
1485  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
1486  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
1487 
1488 
1489  ! Obtaining coefficients of the equation describing the sweep line
1490  a = patch_icpp(patch_id)%normal(1)
1491  b = patch_icpp(patch_id)%normal(2)
1492  c = -a*x_centroid - b*y_centroid
1493 
1494 
1495  ! Initializing the pseudo volume fraction value to 1. The value will
1496  ! be modified as the patch is laid out on the grid, but only in the
1497  ! case that smoothing of the sweep line patch's boundary is enabled.
1498  eta = 1d0
1499 
1500  ! Checking whether the region swept by the line covers a particular
1501  ! cell in the domain and verifying whether the current patch has the
1502  ! permission to write to that cell. If both queries check out, the
1503  ! primitive variables of the current patch are written to this cell.
1504  DO j = 0, n
1505  DO i = 0, m
1506 
1507  IF(patch_icpp(patch_id)%smoothen) THEN
1508  eta = 5d-1 + 5d-1*tanh( smooth_coeff/min(dx,dy) &
1509  * (a*x_cc(i)+b*y_cc(j)+c) &
1510  / sqrt(a**2d0+b**2d0) )
1511  END IF
1512 
1513  IF(( a*x_cc(i) + b*y_cc(j) + c >= 0d0 &
1514  .AND. &
1515  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0))) &
1516  .OR. &
1517  patch_id_fp(i,j,0) == smooth_patch_id ) &
1518  THEN
1519  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1520 
1521  END IF
1522 
1523  END DO
1524  END DO
1525 
1526 
1527  END SUBROUTINE s_sweep_line ! ------------------------------------------
1528 
1529 
1536  SUBROUTINE s_isentropic_vortex(patch_id) ! ----------------------------
1538 
1539  ! Patch identifier
1540  INTEGER, INTENT(IN) :: patch_id
1541 
1542  ! Generic loop iterators
1543  INTEGER :: i,j
1544 
1545 
1546  ! Transferring isentropic vortex patch's centroid and radius info
1547  x_centroid = patch_icpp(patch_id)%x_centroid
1548  y_centroid = patch_icpp(patch_id)%y_centroid
1549  radius = patch_icpp(patch_id)%radius
1550 
1551 
1552  ! Since the isentropic vortex patch does not allow for its boundary
1553  ! to get smoothed, the pseudo volume fraction is set to 1 to ensure
1554  ! that only the current patch contributes to the fluid state in the
1555  ! cells that this patch covers.
1556  eta = 1d0
1557 
1558 
1559  ! Verifying whether the isentropic vortex includes a particular cell
1560  ! and verifying whether the current patch has permission to write to
1561  ! that cell. If both queries work out the primitive variables of the
1562  ! the current patch are assigned to this cell.
1563  DO j = 0, n
1564  DO i = 0, m
1565 
1566  IF( (x_cc(i) - x_centroid)**2d0 &
1567  + (y_cc(j) - y_centroid)**2d0 <= radius**2d0 &
1568  .AND. &
1569  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0)) ) &
1570  THEN
1571 
1572  CALL s_assign_patch_primitive_variables( patch_id, &
1573  i,j,0)
1574 
1575  END IF
1576 
1577  END DO
1578  END DO
1579 
1580 
1581  END SUBROUTINE s_isentropic_vortex ! -----------------------------------
1582 
1586  SUBROUTINE s_1d_analytical(patch_id) ! ---------------------------------
1588  ! Patch identifier
1589  INTEGER, INTENT(IN) :: patch_id
1590 
1591  ! Placeholders for the cell boundary values
1592  REAL(KIND(0d0)) :: a,b,c,d, pi_inf, gamma, lit_gamma
1593 
1594  ! Generic loop iterators
1595  INTEGER :: i,j
1596 
1597  pi_inf = fluid_pp(1)%pi_inf
1598  gamma = fluid_pp(1)%gamma
1599  lit_gamma = (1d0+gamma)/gamma
1600 
1601  ! Transferring the patch's centroid and length information
1602  x_centroid = patch_icpp(patch_id)%x_centroid
1603  length_x = patch_icpp(patch_id)%length_x
1604 
1605  ! Computing the beginning and the end x- and y-coordinates
1606  ! of the patch based on its centroid and lengths
1607  x_boundary%beg = x_centroid - 0.5d0*length_x
1608  x_boundary%end = x_centroid + 0.5d0*length_x
1609 
1610  ! Since the patch doesn't allow for its boundaries to be
1611  ! smoothed out, the pseudo volume fraction is set to 1 to
1612  ! ensure that only the current patch contributes to the fluid
1613  ! state in the cells that this patch covers.
1614  eta = 1d0
1615 
1616  ! Checking whether the line segment covers a particular cell in the
1617  ! domain and verifying whether the current patch has the permission
1618  ! to write to that cell. If both queries check out, the primitive
1619  ! variables of the current patch are assigned to this cell.
1620  DO i = 0, m
1621  IF( x_boundary%beg <= x_cc(i) .AND. &
1622  x_boundary%end >= x_cc(i) .AND. &
1623  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,0,0)) ) THEN
1624 
1625  CALL s_assign_patch_primitive_variables(patch_id, i,0,0)
1626 
1627  !what variables to alter
1628  !bump in pressure
1629  q_prim_vf(e_idx)%sf(i,0,0) = q_prim_vf(e_idx)%sf(i,0,0) * &
1630  ( 1d0 + 0.2d0*dexp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1631 
1632  !bump in void fraction
1633  !q_prim_vf(adv_idx%beg)%sf(i,0,0) = q_prim_vf(adv_idx%beg)%sf(i,0,0) * &
1634  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1635 
1636  !bump in R(x)
1637  !q_prim_vf(adv_idx%end+1)%sf(i,0,0) = q_prim_vf(adv_idx%end+1)%sf(i,0,0) * &
1638  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1639 
1640  !IF (model_eqns == 4) THEN
1641  !reassign density
1642  !IF (num_fluids == 1) THEN
1643  q_prim_vf(1)%sf(i,0,0) = &
1644  (((q_prim_vf(e_idx)%sf(i,0,0) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma)) * &
1645  rhoref*(1d0-q_prim_vf(alf_idx)%sf(i,0,0))
1646  !END IF
1647  !ELSE IF (model_eqns == 2) THEN
1648  !can manually adjust density here
1649  !q_prim_vf(1)%sf(i,0,0) = q_prim_vf(1)%sf(i,0,0) * &
1650  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1651  !END IF
1652  END IF
1653  END DO
1654 
1655 
1656  END SUBROUTINE s_1d_analytical ! ---------------------------------------
1657 
1658 
1659 
1660  SUBROUTINE s_1d_bubble_pulse(patch_id) ! ---------------------------------
1661  ! Description: This patch assigns the primitive variables as analytical
1662  ! functions such that the code can be verified.
1663 
1664  ! Patch identifier
1665  INTEGER, INTENT(IN) :: patch_id
1666 
1667  ! Placeholders for the cell boundary values
1668  REAL(KIND(0d0)) :: fac, a,b,c,d, pi_inf, gamma, lit_gamma
1669 
1670  ! Generic loop iterators
1671  INTEGER :: i,j
1672 
1673  pi_inf = fluid_pp(1)%pi_inf
1674  gamma = fluid_pp(1)%gamma
1675  lit_gamma = (1d0+gamma)/gamma
1676 
1677  ! Transferring the patch's centroid and length information
1678  x_centroid = patch_icpp(patch_id)%x_centroid
1679  length_x = patch_icpp(patch_id)%length_x
1680 
1681 
1682  ! Computing the beginning and the end x- and y-coordinates
1683  ! of the patch based on its centroid and lengths
1684  x_boundary%beg = x_centroid - 0.5d0*length_x
1685  x_boundary%end = x_centroid + 0.5d0*length_x
1686 
1687  ! Since the patch doesn't allow for its boundaries to be
1688  ! smoothed out, the pseudo volume fraction is set to 1 to
1689  ! ensure that only the current patch contributes to the fluid
1690  ! state in the cells that this patch covers.
1691  eta = 1d0
1692 
1693  ! Checking whether the line segment covers a particular cell in the
1694  ! domain and verifying whether the current patch has the permission
1695  ! to write to that cell. If both queries check out, the primitive
1696  ! variables of the current patch are assigned to this cell.
1697  DO i = 0, m
1698  IF( x_boundary%beg <= x_cc(i) .AND. &
1699  x_boundary%end >= x_cc(i) .AND. &
1700  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,0,0)) ) THEN
1701 
1702  CALL s_assign_patch_primitive_variables(patch_id, i,0,0)
1703 
1704  !what variables to alter
1705  !sinusoid in pressure
1706  q_prim_vf(e_idx)%sf(i,0,0) = q_prim_vf(e_idx)%sf(i,0,0) * &
1707  ( 1d0 + 0.1d0*sin(-1d0*(x_cb(i)-x_centroid)*2d0*pi/length_x ) )
1708 
1709  !bump in void fraction
1710  !q_prim_vf(adv_idx%beg)%sf(i,0,0) = q_prim_vf(adv_idx%beg)%sf(i,0,0) * &
1711  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1712 
1713  !bump in R(x)
1714  !q_prim_vf(adv_idx%end+1)%sf(i,0,0) = q_prim_vf(adv_idx%end+1)%sf(i,0,0) * &
1715  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1716 
1717  !IF (model_eqns == 4) THEN
1718  !reassign density
1719  !IF (num_fluids == 1) THEN
1720  q_prim_vf(1)%sf(i,0,0) = &
1721  (((q_prim_vf(e_idx)%sf(i,0,0) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma)) * &
1722  rhoref*(1d0-q_prim_vf(alf_idx)%sf(i,0,0))
1723  !END IF
1724  !ELSE IF (model_eqns == 2) THEN
1725  !can manually adjust density here
1726  !q_prim_vf(1)%sf(i,0,0) = q_prim_vf(1)%sf(i,0,0) * &
1727  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1728  !END IF
1729  END IF
1730  END DO
1731 
1732 
1733  END SUBROUTINE s_1d_bubble_pulse ! ---------------------------------------
1734 
1738  SUBROUTINE s_2d_analytical(patch_id) ! ---------------------------------
1740  INTEGER, INTENT(IN) :: patch_id
1741 
1742  REAL(KIND(0d0)) :: a,b,c,d
1743  REAL(KIND(0d0)) :: pi_inf,gamma,lit_gamma
1744 
1745  INTEGER :: i,j
1746 
1747  pi_inf = fluid_pp(1)%pi_inf
1748  gamma = fluid_pp(1)%gamma
1749  lit_gamma = (1d0+gamma)/gamma
1750 
1751  ! Transferring the patch's centroid and length information
1752  x_centroid = patch_icpp(patch_id)%x_centroid
1753  y_centroid = patch_icpp(patch_id)%y_centroid
1754  length_x = patch_icpp(patch_id)%length_x
1755  length_y = patch_icpp(patch_id)%length_y
1756 
1757  ! Computing the beginning and the end x- and y-coordinates
1758  ! of the patch based on its centroid and lengths
1759  x_boundary%beg = x_centroid - 0.5d0*length_x
1760  x_boundary%end = x_centroid + 0.5d0*length_x
1761  y_boundary%beg = y_centroid - 0.5d0*length_y
1762  y_boundary%end = y_centroid + 0.5d0*length_y
1763 
1764  ! Since the patch doesn't allow for its boundaries to be
1765  ! smoothed out, the pseudo volume fraction is set to 1 to
1766  ! ensure that only the current patch contributes to the fluid
1767  ! state in the cells that this patch covers.
1768  eta = 1d0
1769 
1770  ! Checking whether the patch covers a particular cell in the
1771  ! domain and verifying whether the current patch has the
1772  ! permission to write to that cell. If both queries check out,
1773  ! the primitive variables of the current patch are assigned
1774  ! to this cell.
1775  DO j = 0, n
1776  DO i = 0, m
1777  IF ( x_boundary%beg <= x_cc(i) .AND. &
1778  x_boundary%end >= x_cc(i) .AND. &
1779  y_boundary%beg <= y_cc(j) .AND. &
1780  y_boundary%end >= y_cc(j) .AND. &
1781  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,0))) THEN
1782 
1783  CALL s_assign_patch_primitive_variables(patch_id, i,j,0)
1784 
1785  !what variables to alter
1786  !x-y bump in pressure
1787  !q_prim_vf(E_idx)%sf(i,j,0) = q_prim_vf(E_idx)%sf(i,j,0) * &
1788  ! ( 1d0 + 0.2d0*dexp(-1d0*((x_cb(i)-x_centroid)**2.d0 + (y_cb(j)-y_centroid)**2.d0)/(2.d0*0.005d0)) )
1789 
1790  !x-bump
1791  q_prim_vf(e_idx)%sf(i,j,0) = q_prim_vf(e_idx)%sf(i,j,0) * &
1792  ( 1d0 + 0.2d0*dexp(-1d0*((x_cb(i)-x_centroid)**2.d0)/(2.d0*0.005d0)) )
1793 
1794  !bump in void fraction
1795  !q_prim_vf(adv_idx%beg)%sf(i,j,0) = q_prim_vf(adv_idx%beg)%sf(i,j,0) * &
1796  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0 + (y_cb(j)-y_centroid)**2.d0)/(2.d0*0.005d0)) )
1797 
1798  !bump in R(x)
1799  !q_prim_vf(adv_idx%end+1)%sf(i,j,0) = q_prim_vf(adv_idx%end+1)%sf(i,j,0) * &
1800  ! ( 1d0 + 0.2d0*exp(-1d0*((x_cb(i)-x_centroid)**2.d0 + (y_cb(j)-y_centroid)**2.d0)/(2.d0*0.005d0)) )
1801 
1802  !reassign density
1803  q_prim_vf(1)%sf(i,j,0) = &
1804  (((q_prim_vf(e_idx)%sf(i,j,0) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma)) * &
1805  rhoref*(1d0-q_prim_vf(alf_idx)%sf(i,j,0))
1806 
1807  ! ================================================================================
1808 
1809  ! Sinusoidal initial condition for all flow variables =============================
1810 
1811  ! Cell-center values
1812 ! a = 0d0
1813 ! b = 0d0
1814 ! c = 0d0
1815 ! d = 0d0
1816 ! q_prim_vf(adv_idx%beg)%sf(i,j,0) = SIN(x_cc(i)) * SIN(y_cc(j))
1817 ! q_prim_vf(1)%sf(i,j,0) = q_prim_vf(adv_idx%beg)%sf(i,j,0) * 1d0
1818 ! q_prim_vf(cont_idx%end)%sf(i,j,0) = (1d0 - q_prim_vf(adv_idx%beg)%sf(i,j,0)) * 1d0
1819 ! q_prim_vf(mom_idx%beg)%sf(i,j,0) = SIN(x_cc(i))
1820 ! q_prim_vf(mom_idx%end)%sf(i,j,0) = SIN(y_cc(j))
1821 ! q_prim_vf(E_idx)%sf(i,j,0) = 1d0
1822 
1823  ! Cell-average values
1824 ! a = x_cc(i) - 5d-1*dx ! x-beg
1825 ! b = x_cc(i) + 5d-1*dx ! x-end
1826 ! c = y_cc(j) - 5d-1*dy ! y-beg
1827 ! d = y_cc(j) + 5d-1*dy ! y-end
1828 ! q_prim_vf(adv_idx%beg)%sf(i,j,0) = 1d0/((b-a)*(d-c)) * &
1829 ! (COS(a)*COS(c) - COS(a)*COS(d) - COS(b)*COS(c) + COS(b)*COS(d))
1830 ! q_prim_vf(1)%sf(i,j,0) = q_prim_vf(adv_idx%beg)%sf(i,j,0) * 1d0
1831 ! q_prim_vf(cont_idx%end)%sf(i,j,0) = (1d0 - q_prim_vf(adv_idx%beg)%sf(i,j,0)) * 1d0
1832 ! q_prim_vf(mom_idx%beg)%sf(i,j,0) = (COS(a) - COS(b))/(b-a)
1833 ! q_prim_vf(mom_idx%end)%sf(i,j,0) = (COS(c) - COS(d))/(d-c)
1834 ! q_prim_vf(E_idx)%sf(i,j,0) = 1d0
1835  ! ================================================================================
1836 
1837  ! Initial pressure profile smearing for bubble collapse case of Tiwari (2013) ====
1838  !IF(( (x_cc(i))**2d0 &
1839  ! + (y_cc(j))**2d0 <= 1d0**2d0)) THEN
1840  ! q_prim_vf(E_idx)%sf(i,j,0) = 1d5 / 25d0
1841  !ELSE
1842  ! q_prim_vf(E_idx)%sf(i,j,0) = 1d5 + 1d0/SQRT(x_cc(i)**2d0+y_cc(j)**2d0) &
1843  ! * ((1d5/25d0) - 1d5)
1844  !END IF
1845  ! ================================================================================
1846 
1847  END IF
1848  END DO
1849  END DO
1850 
1851  END SUBROUTINE s_2d_analytical ! ---------------------------------------
1852 
1853 
1854 
1858  SUBROUTINE s_3d_analytical(patch_id) ! ---------------------------------
1860 
1861  INTEGER, INTENT(IN) :: patch_id
1862  REAL(KIND(0d0)) :: pi_inf, gamma, lit_gamma
1863 
1864  INTEGER :: i,j,k
1865 
1866  pi_inf = fluid_pp(1)%pi_inf
1867  gamma = fluid_pp(1)%gamma
1868  lit_gamma = (1d0+gamma)/gamma
1869 
1870  ! Transferring the patch's centroid and length information
1871  x_centroid = patch_icpp(patch_id)%x_centroid
1872  y_centroid = patch_icpp(patch_id)%y_centroid
1873  z_centroid = patch_icpp(patch_id)%z_centroid
1874  length_x = patch_icpp(patch_id)%length_x
1875  length_y = patch_icpp(patch_id)%length_y
1876  length_z = patch_icpp(patch_id)%length_z
1877 
1878  ! Computing the beginning and the end x-, y- and z-coordinates of
1879  ! the patch based on its centroid and lengths
1880  x_boundary%beg = x_centroid - 0.5d0*length_x
1881  x_boundary%end = x_centroid + 0.5d0*length_x
1882  y_boundary%beg = y_centroid - 0.5d0*length_y
1883  y_boundary%end = y_centroid + 0.5d0*length_y
1884  z_boundary%beg = z_centroid - 0.5d0*length_z
1885  z_boundary%end = z_centroid + 0.5d0*length_z
1886 
1887  ! Since the analytical patch does not allow for its boundaries to get
1888  ! smoothed out, the pseudo volume fraction is set to 1 to make sure
1889  ! that only the current patch contributes to the fluid state in the
1890  ! cells that this patch covers.
1891  eta = 1d0
1892 
1893  ! Checking whether the patch covers a particular cell in the domain
1894  ! and verifying whether the current patch has permission to write to
1895  ! to that cell. If both queries check out, the primitive variables
1896  ! of the current patch are assigned to this cell.
1897  DO k = 0, p
1898  DO j = 0, n
1899  DO i = 0, m
1900 
1901  IF (grid_geometry == 3) THEN
1902  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
1903  ELSE
1904  cart_y = y_cc(j)
1905  cart_z = z_cc(k)
1906  END IF
1907 
1908  IF( x_boundary%beg <= x_cc(i) .AND. &
1909  x_boundary%end >= x_cc(i) .AND. &
1910  y_boundary%beg <= cart_y .AND. &
1911  y_boundary%end >= cart_y .AND. &
1912  z_boundary%beg <= cart_z .AND. &
1913  z_boundary%end >= cart_z &
1914  .AND. &
1915  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k))) &
1916  THEN
1917 
1918  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
1919 
1920  !gaussian ball
1921  !what variables to alter
1922  !bump in pressure
1923  q_prim_vf(e_idx)%sf(i,j,k) = q_prim_vf(e_idx)%sf(i,j,k) * &
1924  ( 1d0 + 0.2d0*exp(-1d0 * &
1925  ((x_cb(i)-x_centroid)**2.d0 + (y_cb(j)-y_centroid)**2.d0 + (z_cb(k)-z_centroid)**2.d0) &
1926  /(2.d0*0.005d0)) )
1927 
1928  !bump in void fraction
1929  q_prim_vf(adv_idx%beg)%sf(i,j,k) = q_prim_vf(adv_idx%beg)%sf(i,j,k) * &
1930  ( 1d0 + 0.2d0*exp(-1d0 * &
1931  ((x_cb(i)-x_centroid)**2.d0 + (y_cb(j)-y_centroid)**2.d0 + (z_cb(k)-z_centroid)**2.d0) &
1932  /(2.d0*0.005d0)) )
1933 
1934  !bump in R(x)
1935  q_prim_vf(adv_idx%end+1)%sf(i,j,k) = q_prim_vf(adv_idx%end+1)%sf(i,j,k) * &
1936  ( 1d0 + 0.2d0*exp(-1d0 * &
1937  ((x_cb(i)-x_centroid)**2.d0 + (y_cb(j)-y_centroid)**2.d0 + (z_cb(k)-z_centroid)**2.d0) &
1938  /(2.d0*0.005d0)) )
1939 
1940  !reassign density
1941  q_prim_vf(1)%sf(i,j,k) = &
1942  (((q_prim_vf(e_idx)%sf(i,j,k) + pi_inf)/(pref + pi_inf))**(1d0/lit_gamma)) * &
1943  rhoref*(1d0-q_prim_vf(e_idx+1)%sf(i,j,k))
1944 
1945  ! ================================================================================
1946 
1947  ! Constant x-velocity in cylindrical grid ========================================
1948 ! q_prim_vf(cont_idx%beg )%sf(i,j,k) = 1d0
1949 ! q_prim_vf(cont_idx%end )%sf(i,j,k) = 0d0
1950 ! q_prim_vf(mom_idx%beg )%sf(i,j,k) = 0d0
1951 ! q_prim_vf(mom_idx%beg+1)%sf(i,j,k) = COS(z_cc(k))
1952 ! q_prim_vf(mom_idx%end )%sf(i,j,k) = -SIN(z_cc(k))
1953 ! q_prim_vf(E_idx )%sf(i,j,k) = 1d0
1954 ! q_prim_vf(adv_idx%beg )%sf(i,j,k) = 1d0
1955  ! ================================================================================
1956 
1957  ! Couette flow in cylindrical grid ===============================================
1958  !q_prim_vf(cont_idx%beg )%sf(i,j,k) = 1d0
1959  !q_prim_vf(cont_idx%end )%sf(i,j,k) = 0d0
1960  !q_prim_vf(mom_idx%beg )%sf(i,j,k) = 0d0
1961  !q_prim_vf(mom_idx%beg+1)%sf(i,j,k) = y_cc(j)*COS(z_cc(k))*SIN(z_cc(k))
1962  !q_prim_vf(mom_idx%end )%sf(i,j,k) = -y_cc(j)*SIN(z_cc(k))**2d0
1963  !q_prim_vf(E_idx )%sf(i,j,k) = 1d0
1964  !q_prim_vf(adv_idx%beg )%sf(i,j,k) = 1d0
1965  ! ================================================================================
1966 
1967  END IF
1968 
1969  END DO
1970  END DO
1971  END DO
1972 
1973  END SUBROUTINE s_3d_analytical ! ---------------------------------------
1974 
1975 
1976 
1980  SUBROUTINE s_spherical_harmonic(patch_id) ! ----------------------------
1982  INTEGER, INTENT(IN) :: patch_id
1983 
1984  INTEGER :: i,j,k
1985 
1986  COMPLEX(KIND(0d0)) :: cmplx_i = (0d0,1d0)
1987  COMPLEX(KIND(0d0)) :: H
1988 
1989  ! Transferring the patch's centroid and radius information
1990  x_centroid = patch_icpp(patch_id)%x_centroid
1991  y_centroid = patch_icpp(patch_id)%y_centroid
1992  z_centroid = patch_icpp(patch_id)%z_centroid
1993  radius = patch_icpp(patch_id)%radius
1994  epsilon = patch_icpp(patch_id)%epsilon
1995  beta = patch_icpp(patch_id)%beta
1996 
1997  ! Since the analytical patch does not allow for its boundaries to get
1998  ! smoothed out, the pseudo volume fraction is set to 1 to make sure
1999  ! that only the current patch contributes to the fluid state in the
2000  ! cells that this patch covers.
2001  eta = 1d0
2002 
2003  ! Checking whether the patch covers a particular cell in the domain
2004  ! and verifying whether the current patch has permission to write to
2005  ! to that cell. If both queries check out, the primitive variables
2006  ! of the current patch are assigned to this cell.
2007  DO k = 0, p
2008  DO j = 0, n
2009  DO i = 0, m
2010 
2011  IF (grid_geometry == 3) THEN
2012  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
2013  ELSE
2014  cart_y = y_cc(j)
2015  cart_z = z_cc(k)
2016  END IF
2017 
2018  IF(( ( x_cc(i) - x_centroid)**2d0 &
2019  + ( cart_y - y_centroid)**2d0 &
2020  + ( cart_z - z_centroid)**2d0 <= radius**2d0 &
2021  .AND. &
2022  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k)))) &
2023  THEN
2024 
2025  CALL s_convert_cylindrical_to_spherical_coord(x_cc(i),y_cc(j))
2026 
2027  IF (epsilon == 1d0) THEN
2028  IF (beta == 0d0) THEN
2029  h = 5d-1*sqrt(3d0/pi)*cos(sph_phi)
2030  ELSEIF (beta == 1d0) THEN
2031  h = -5d-1*sqrt(3d0/(2d0*pi))*exp(cmplx_i*z_cc(k))*sin(sph_phi)
2032  END IF
2033  ELSEIF (epsilon == 2d0) THEN
2034  IF (beta == 0d0) THEN
2035  h = 25d-2*sqrt(5d0/pi)*(3d0*cos(sph_phi)**2d0 - 1d0)
2036  ELSEIF (beta == 1d0) THEN
2037  h = -5d-1*sqrt(15d0/(2d0*pi))*exp(cmplx_i*z_cc(k))*sin(sph_phi)*cos(sph_phi)
2038  ELSEIF (beta == 2d0) THEN
2039  h = 25d-2*sqrt(15d0/(2d0*pi))*exp(2d0*cmplx_i*z_cc(k))*sin(sph_phi)**2d0
2040  END IF
2041  ELSEIF (epsilon == 3d0) THEN
2042  IF (beta == 0d0) THEN
2043  h = 25d-2*sqrt(7d0/pi)*(5d0*cos(sph_phi)**3d0 - 3d0*cos(sph_phi))
2044  ELSEIF (beta == 1d0) THEN
2045  h = -125d-3*sqrt(21d0/pi)*exp(cmplx_i*z_cc(k))*sin(sph_phi) * &
2046  (5d0*cos(sph_phi)**2d0 - 1d0)
2047  ELSEIF (beta == 2d0) THEN
2048  h = 25d-2*sqrt(105d0/(2d0*pi))*exp(2d0*cmplx_i*z_cc(k)) * &
2049  sin(sph_phi)**2d0*cos(sph_phi)
2050  ELSEIF (beta == 3d0) THEN
2051  h = -125d-3*sqrt(35d0/pi)*exp(3d0*cmplx_i*z_cc(k))*sin(sph_phi)**3d0
2052  END IF
2053  ELSEIF (epsilon == 4d0) THEN
2054  IF (beta == 0d0) THEN
2055  h = 3d0/16d0*sqrt(1d0/pi)*(35d0*cos(sph_phi)**4d0 - &
2056  3d1*cos(sph_phi)**2d0 + 3d0)
2057  ELSEIF (beta == 1d0) THEN
2058  h = -3d0/8d0*sqrt(5d0/pi)*exp(cmplx_i*z_cc(k))* &
2059  sin(sph_phi)*(7d0*cos(sph_phi)**3d0 - 3d0*cos(sph_phi))
2060  ELSEIF (beta == 2d0) THEN
2061  h = 3d0/8d0*sqrt(5d0/(2d0*pi))*exp(2d0*cmplx_i*z_cc(k))* &
2062  sin(sph_phi)**2d0*(7d0*cos(sph_phi)**2d0 - 1d0)
2063  ELSEIF (beta == 3d0) THEN
2064  h = -3d0/8d0*sqrt(35d0/pi)*exp(3d0*cmplx_i*z_cc(k))* &
2065  sin(sph_phi)**3d0*cos(sph_phi)
2066  ELSEIF (beta == 4d0) THEN
2067  h = 3d0/16d0*sqrt(35d0/(2d0*pi))*exp(4d0*cmplx_i*z_cc(k))* &
2068  sin(sph_phi)**4d0
2069  END IF
2070  ELSEIF (epsilon == 5d0) THEN
2071  IF (beta == 0d0) THEN
2072  h = 1d0/16d0*sqrt(11d0/pi)*(63d0*cos(sph_phi)**5d0 - &
2073  7d1*cos(sph_phi)**3d0 + 15d0*cos(sph_phi))
2074  ELSEIF (beta == 1d0) THEN
2075  h = -1d0/16d0*sqrt(165d0/(2d0*pi))*exp(cmplx_i*z_cc(k))* &
2076  sin(sph_phi)*(21d0*cos(sph_phi)**4d0 - 14d0*cos(sph_phi)**2d0 + 1d0)
2077  ELSEIF (beta == 2d0) THEN
2078  h = 125d-3*sqrt(1155d0/(2d0*pi))*exp(2d0*cmplx_i*z_cc(k)) * &
2079  sin(sph_phi)**2d0*(3d0*cos(sph_phi)**3d0 - cos(sph_phi))
2080  ELSEIF (beta == 3d0) THEN
2081  h = -1d0/32d0*sqrt(385d0/pi)*exp(3d0*cmplx_i*z_cc(k) )* &
2082  sin(sph_phi)**3d0*(9d0*cos(sph_phi)**2d0 - 1d0)
2083  ELSEIF (beta == 4d0) THEN
2084  h = 3d0/16d0*sqrt(385d0/(2d0*pi))*exp(4d0*cmplx_i*z_cc(k)) * &
2085  sin(sph_phi)**4d0*cos(sph_phi)
2086  ELSEIF (beta == 5d0) THEN
2087  h = -3d0/32d0*sqrt(77d0/pi)*exp(5d0*cmplx_i*z_cc(k)) * &
2088  sin(sph_phi)**5d0
2089  END IF
2090  END IF
2091 
2092  q_prim_vf(adv_idx%beg)%sf(i,j,k) = 1d0-abs(REAL(h,kind(0d0)))
2093 
2094  END IF
2095 
2096  END DO
2097  END DO
2098  END DO
2099 
2100  END SUBROUTINE s_spherical_harmonic ! ----------------------------------
2101 
2102 
2109  SUBROUTINE s_sphere(patch_id) ! ----------------------------------------
2111  INTEGER, INTENT(IN) :: patch_id
2112 
2113  ! Generic loop iterators
2114  INTEGER :: i,j,k
2115 
2116 
2117  REAL(KIND(0d0)) :: radius_pressure, pressure_bubble, pressure_inf
2120 
2121 
2122  ! Transferring spherical patch's radius, centroid, smoothing patch
2123  ! identity and smoothing coefficient information
2124  x_centroid = patch_icpp(patch_id)%x_centroid
2125  y_centroid = patch_icpp(patch_id)%y_centroid
2126  z_centroid = patch_icpp(patch_id)%z_centroid
2127  radius = patch_icpp(patch_id)%radius
2128  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
2129  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
2130 
2131 
2132  ! Initializing the pseudo volume fraction value to 1. The value will
2133  ! be modified as the patch is laid out on the grid, but only in the
2134  ! case that smoothing of the spherical patch's boundary is enabled.
2135  eta = 1d0
2136 
2137 
2138  ! Checking whether the sphere covers a particular cell in the domain
2139  ! and verifying whether the current patch has permission to write to
2140  ! that cell. If both queries check out, the primitive variables of
2141  ! the current patch are assigned to this cell.
2142  DO k = 0, p
2143  DO j = 0, n
2144  DO i = 0, m
2145 
2146  IF (grid_geometry == 3) THEN
2147  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
2148  ELSE
2149  cart_y = y_cc(j)
2150  cart_z = z_cc(k)
2151  END IF
2152 
2153  IF(patch_icpp(patch_id)%smoothen) THEN
2154 
2155  eta = tanh( smooth_coeff / min(dx,dy,dz) * &
2156  ( sqrt( ( x_cc(i) - x_centroid)**2d0 &
2157  + ( cart_y - y_centroid)**2d0 &
2158  + ( cart_z - z_centroid)**2d0 ) &
2159  - radius ) ) * (-0.5d0) + 0.5d0
2160 
2161  END IF
2162 
2163  IF(( ( x_cc(i) - x_centroid)**2d0 &
2164  + ( cart_y - y_centroid)**2d0 &
2165  + ( cart_z - z_centroid)**2d0 <= radius**2d0 &
2166  .AND. &
2167  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k))) &
2168  .OR. &
2169  patch_id_fp(i,j,k)== smooth_patch_id ) &
2170  THEN
2171 
2172  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
2173 
2174  END IF
2175 
2176  ! Initialization of the pressure field that corresponds to the bubble-collapse
2177  !! test case found in Tiwari et al. (2013)
2178  ! radius_pressure = SQRT(x_cc(i)**2d0) ! 1D
2179  ! radius_pressure = SQRT(x_cc(i)**2d0 + cart_y**2d0) ! 2D
2180  ! radius_pressure = SQRT(x_cc(i)**2d0 + cart_y**2d0 + cart_z**2d0) ! 3D
2181  ! pressure_bubble = 1.E+04
2182  ! pressure_inf = 1.E+05
2183  ! q_prim_vf(E_idx)%sf(i,j,k) = pressure_inf + radius / radius_pressure * (pressure_bubble - pressure_inf)
2184  !
2185  ! IF( (( x_cc(i) - x_centroid)**2d0 &
2186  ! + ( cart_y - y_centroid)**2d0 &
2187  ! + ( cart_z - z_centroid)**2d0) <= radius**2d0) &
2188  ! THEN
2189  !
2190  ! q_prim_vf(E_idx)%sf(i,j,k) = pressure_bubble
2191  !
2192  ! END IF
2193 
2194  END DO
2195  END DO
2196  END DO
2197 
2198 
2199  END SUBROUTINE s_sphere ! ----------------------------------------------
2200 
2201 
2211  SUBROUTINE s_cuboid(patch_id) ! ----------------------------------------
2213  INTEGER, INTENT(IN) :: patch_id
2214 
2215  INTEGER :: i,j,k
2216 
2217 
2218  ! Transferring the cuboid's centroid and length information
2219  x_centroid = patch_icpp(patch_id)%x_centroid
2220  y_centroid = patch_icpp(patch_id)%y_centroid
2221  z_centroid = patch_icpp(patch_id)%z_centroid
2222  length_x = patch_icpp(patch_id)%length_x
2223  length_y = patch_icpp(patch_id)%length_y
2224  length_z = patch_icpp(patch_id)%length_z
2225 
2226 
2227  ! Computing the beginning and the end x-, y- and z-coordinates of
2228  ! the cuboid based on its centroid and lengths
2229  x_boundary%beg = x_centroid - 0.5d0*length_x
2230  x_boundary%end = x_centroid + 0.5d0*length_x
2231  y_boundary%beg = y_centroid - 0.5d0*length_y
2232  y_boundary%end = y_centroid + 0.5d0*length_y
2233  z_boundary%beg = z_centroid - 0.5d0*length_z
2234  z_boundary%end = z_centroid + 0.5d0*length_z
2235 
2236 
2237  ! Since the cuboidal patch does not allow for its boundaries to get
2238  ! smoothed out, the pseudo volume fraction is set to 1 to make sure
2239  ! that only the current patch contributes to the fluid state in the
2240  ! cells that this patch covers.
2241  eta = 1d0
2242 
2243 
2244  ! Checking whether the cuboid covers a particular cell in the domain
2245  ! and verifying whether the current patch has permission to write to
2246  ! to that cell. If both queries check out, the primitive variables
2247  ! of the current patch are assigned to this cell.
2248  DO k = 0, p
2249  DO j = 0, n
2250  DO i = 0, m
2251 
2252  IF (grid_geometry == 3) THEN
2253  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
2254  ELSE
2255  cart_y = y_cc(j)
2256  cart_z = z_cc(k)
2257  END IF
2258 
2259  IF( x_boundary%beg <= x_cc(i) .AND. &
2260  x_boundary%end >= x_cc(i) .AND. &
2261  y_boundary%beg <= cart_y .AND. &
2262  y_boundary%end >= cart_y .AND. &
2263  z_boundary%beg <= cart_z .AND. &
2264  z_boundary%end >= cart_z &
2265  .AND. &
2266  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k))) &
2267  THEN
2268 
2269  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
2270 
2271  END IF
2272  END DO
2273  END DO
2274  END DO
2275 
2276 
2277  END SUBROUTINE s_cuboid ! ----------------------------------------------
2278 
2279 
2289  SUBROUTINE s_cylinder(patch_id) ! --------------------------------------
2291  INTEGER, INTENT(IN) :: patch_id
2292 
2293  INTEGER :: i,j,k
2294 
2295 
2296  ! Transferring the cylindrical patch's centroid, length, radius,
2297  ! smoothing patch identity and smoothing coefficient information
2298  x_centroid = patch_icpp(patch_id)%x_centroid
2299  y_centroid = patch_icpp(patch_id)%y_centroid
2300  z_centroid = patch_icpp(patch_id)%z_centroid
2301  length_x = patch_icpp(patch_id)%length_x
2302  length_y = patch_icpp(patch_id)%length_y
2303  length_z = patch_icpp(patch_id)%length_z
2304  radius = patch_icpp(patch_id)%radius
2305  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
2306  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
2307 
2308 
2309  ! Computing the beginning and the end x-, y- and z-coordinates of
2310  ! the cylinder based on its centroid and lengths
2311  x_boundary%beg = x_centroid - 0.5d0*length_x
2312  x_boundary%end = x_centroid + 0.5d0*length_x
2313  y_boundary%beg = y_centroid - 0.5d0*length_y
2314  y_boundary%end = y_centroid + 0.5d0*length_y
2315  z_boundary%beg = z_centroid - 0.5d0*length_z
2316  z_boundary%end = z_centroid + 0.5d0*length_z
2317 
2318 
2319  ! Initializing the pseudo volume fraction value to 1. The value will
2320  ! be modified as the patch is laid out on the grid, but only in the
2321  ! case that smearing of the cylindrical patch's boundary is enabled.
2322  eta = 1d0
2323 
2324 
2325  ! Checking whether the cylinder covers a particular cell in the
2326  ! domain and verifying whether the current patch has the permission
2327  ! to write to that cell. If both queries check out, the primitive
2328  ! variables of the current patch are assigned to this cell.
2329  DO k = 0, p
2330  DO j = 0, n
2331  DO i = 0, m
2332 
2333  IF (grid_geometry == 3) THEN
2334  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
2335  ELSE
2336  cart_y = y_cc(j)
2337  cart_z = z_cc(k)
2338  END IF
2339 
2340  IF(patch_icpp(patch_id)%smoothen) THEN
2341 
2342  IF(length_x /= dflt_real) THEN
2343  eta = tanh( smooth_coeff / min(dy,dz) * &
2344  ( sqrt( ( cart_y - y_centroid)**2d0 &
2345  + ( cart_z - z_centroid)**2d0 ) &
2346  - radius ) ) * (-0.5d0) + 0.5d0
2347  ELSEIF(length_y /= dflt_real) THEN
2348  eta = tanh( smooth_coeff / min(dx,dz) * &
2349  ( sqrt( (x_cc(i) - x_centroid)**2d0 &
2350  + ( cart_z - z_centroid)**2d0 ) &
2351  - radius ) ) * (-0.5d0) + 0.5d0
2352  ELSE
2353  eta = tanh( smooth_coeff / min(dx,dy) * &
2354  ( sqrt( (x_cc(i) - x_centroid)**2d0 &
2355  + ( cart_y - y_centroid)**2d0 ) &
2356  - radius ) ) * (-0.5d0) + 0.5d0
2357  END IF
2358 
2359  END IF
2360 
2361  IF(((( length_x /= dflt_real .AND. &
2362  ( cart_y - y_centroid)**2d0 &
2363  + ( cart_z - z_centroid)**2d0 <= radius**2d0 .AND. &
2364  x_boundary%beg <= x_cc(i) .AND. &
2365  x_boundary%end >= x_cc(i) ) &
2366  .OR. &
2367  ( length_y /= dflt_real .AND. &
2368  (x_cc(i) - x_centroid)**2d0 &
2369  + (cart_z - z_centroid)**2d0 <= radius**2d0 .AND. &
2370  y_boundary%beg <= cart_y .AND. &
2371  y_boundary%end >= cart_y ) &
2372  .OR. &
2373  ( length_z /= dflt_real .AND. &
2374  (x_cc(i) - x_centroid)**2d0 &
2375  + ( cart_y - y_centroid)**2d0 <= radius**2d0 .AND. &
2376  z_boundary%beg <= cart_z .AND. &
2377  z_boundary%end >= cart_z )) &
2378  .AND. &
2379  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k)))&
2380  .OR. &
2381  patch_id_fp(i,j,k)== smooth_patch_id )&
2382  THEN
2383 
2384  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
2385 
2386  END IF
2387 
2388  END DO
2389  END DO
2390  END DO
2391 
2392 
2393  END SUBROUTINE s_cylinder ! --------------------------------------------
2394 
2395 
2404  SUBROUTINE s_sweep_plane(patch_id) ! -----------------------------------
2406  INTEGER, INTENT(IN) :: patch_id
2407 
2408  INTEGER :: i,j,k
2409 
2410  ! Transferring the centroid information of the plane to be swept
2411  x_centroid = patch_icpp(patch_id)%x_centroid
2412  y_centroid = patch_icpp(patch_id)%y_centroid
2413  z_centroid = patch_icpp(patch_id)%z_centroid
2414  smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
2415  smooth_coeff = patch_icpp(patch_id)%smooth_coeff
2416 
2417 
2418  ! Obtaining coefficients of the equation describing the sweep plane
2419  a = patch_icpp(patch_id)%normal(1)
2420  b = patch_icpp(patch_id)%normal(2)
2421  c = patch_icpp(patch_id)%normal(3)
2423 
2424 
2425  ! Initializing the pseudo volume fraction value to 1. The value will
2426  ! be modified as the patch is laid out on the grid, but only in the
2427  ! case that smearing of the sweep plane patch's boundary is enabled.
2428  eta = 1d0
2429 
2430 
2431  ! Checking whether the region swept by the plane covers a particular
2432  ! cell in the domain and verifying whether the current patch has the
2433  ! permission to write to that cell. If both queries check out, the
2434  ! primitive variables of the current patch are written to this cell.
2435  DO k = 0, p
2436  DO j = 0, n
2437  DO i = 0, m
2438 
2439  IF (grid_geometry == 3) THEN
2440  CALL s_convert_cylindrical_to_cartesian_coord(y_cc(j),z_cc(k))
2441  ELSE
2442  cart_y = y_cc(j)
2443  cart_z = z_cc(k)
2444  END IF
2445 
2446  IF(patch_icpp(patch_id)%smoothen) THEN
2447  eta = 5d-1 + 5d-1*tanh( smooth_coeff/min(dx,dy,dz) &
2448  * ( a*x_cc(i) + &
2449  b*cart_y + &
2450  c*cart_z + d ) &
2451  / sqrt(a**2d0+b**2d0+c**2d0) )
2452  END IF
2453 
2454  IF(( a*x_cc(i) + b*cart_y + c*cart_z + d >= 0d0 &
2455  .AND. &
2456  patch_icpp(patch_id)%alter_patch(patch_id_fp(i,j,k))) &
2457  .OR. &
2458  patch_id_fp(i,j,k)== smooth_patch_id ) &
2459  THEN
2460 
2461  CALL s_assign_patch_primitive_variables(patch_id, i,j,k)
2462 
2463  END IF
2464 
2465  END DO
2466  END DO
2467  END DO
2468 
2469 
2470  END SUBROUTINE s_sweep_plane ! -----------------------------------------
2471 
2472 
2473 
2474 
2476  SUBROUTINE s_finalize_initial_condition_module() ! ---------------------
2478  INTEGER :: i
2479 
2480 
2481  ! Dellocating the primitive and conservative variables
2482  DO i = 1, sys_size
2483  DEALLOCATE(q_prim_vf(i)%sf)
2484  DEALLOCATE(q_cons_vf(i)%sf)
2485  END DO
2486 
2487  DEALLOCATE(q_prim_vf)
2488  DEALLOCATE(q_cons_vf)
2489 
2490 
2491  ! Deallocating the patch identities bookkeeping variable
2492  DEALLOCATE(patch_id_fp)
2493 
2494 
2495  ! Nullifying procedure pointer to the subroutine assigning either
2496  ! the patch mixture or species primitive variables to a cell in the
2497  ! computational domain
2499 
2500 
2501  END SUBROUTINE s_finalize_initial_condition_module ! -------------------
2502 
2503 
2504 END MODULE m_initial_condition
real(kind(0d0)), dimension(num_fluids_max) fluid_rho
real(kind(0d0)), dimension(:), allocatable v0
subroutine s_varcircle(patch_id)
The varcircle patch is a 2D geometry that may be used . It generatres an annulus. ...
integer, dimension(:,:,:), allocatable patch_id_fp
Bookkepping variable used to track the patch identities (id) associated with each of the cells in the...
subroutine s_sphere(patch_id)
The spherical patch is a 3D geometry that may be used, for example, in creating a bubble or a droplet...
real(kind(0d0)) x_centroid
subroutine s_circle(patch_id)
The circular patch is a 2D geometry that may be used, for example, in creating a bubble or a droplet...
subroutine s_1d_bubble_pulse(patch_id)
subroutine s_cuboid(patch_id)
The cuboidal patch is a 3D geometry that may be used, for example, in creating a solid boundary...
integer perturb_sph_fluid
Fluid to be perturbed with perturb_sph flag.
type(physical_parameters), dimension(num_fluids_max) fluid_pp
Database of the physical parameters of each of the fluids that is present in the flow. These include the stiffened gas equation of state parameters, the Reynolds numbers and the Weber numbers.
integer pi_inf_idx
Index of liquid stiffness func. eqn.
subroutine s_1d_analytical(patch_id)
This patch assigns the primitive variables as analytical functions such that the code can be verified...
real(kind(0d0)), dimension(:), allocatable y_cc
This module contains all of the parameters characterizing the computational domain, simulation algorithm, initial condition and the stiffened equation of state.
real(kind(0d0)) function f_r(myth, offset, a)
Archimedes spiral function.
Derived type annexing a scalar field (SF)
real(kind(0d0)) d
When a line or a plane sweep patch geometry is employed, these variables represent the coefficients a...
subroutine s_convert_species_to_mixture_variables(q_vf, j, k, l, rho, gamma, pi_inf)
This subroutine is designed for the volume fraction model and provided a set of either conservative o...
subroutine s_2d_analytical(patch_id)
This patch assigns the primitive variables as analytical functions such that the code can be verified...
subroutine s_spiral(patch_id)
The spiral patch is a 2D geometry that may be used, The geometry of the patch is well-defined when it...
Abstract interface to the two subroutines that assign the patch primitive variables, either mixture or species, depending on the subroutine, to a particular cell in the computational domain.
type(int_bounds_info) cont_idx
Indexes of first & last continuity eqns.
subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
type(scalar_field), dimension(:), allocatable q_cons_vf
conservative variables
subroutine s_ellipse(patch_id)
The elliptical patch is a 2D geometry. The geometry of the patch is well-defined when its centroid an...
This file contains the definitions of all of the custom-defined types used in the pre-process code...
subroutine s_convert_primitive_to_conservative_variables(q_prim_vf, q_cons_vf)
Converts the primitive variables to the conservative ones. Used when initializing patches...
type(int_bounds_info) adv_idx
Indexes of first & last advection eqns.
integer p
Number of cells in the x-, y- and z-coordinate directions.
real(kind(0d0)), parameter dflt_real
Default real value.
subroutine s_isentropic_vortex(patch_id)
The isentropic vortex is a 2D geometry that may be used, for example, to generate an isentropic flow ...
integer alf_idx
Index of void fraction.
integer num_fluids
Number of different fluids present in the flow.
real(kind(0d0)), dimension(:), allocatable mass_v0
subroutine s_convert_cylindrical_to_spherical_coord(cyl_x, cyl_y)
real(kind(0d0)), dimension(:), allocatable x_cc
logical old_ic
Use existing IC data.
integer sys_size
Number of unknowns in the system of equations.
logical mpp_lim
Alpha limiter.
real(kind(0d0)) smooth_coeff
These variables are analogous in both meaning and use to the similarly named components in the ic_pat...
subroutine s_rectangle(patch_id)
The rectangular patch is a 2D geometry that may be used, for example, in creating a solid boundary...
integer model_eqns
Multicomponent flow model.
type(bub_bounds_info) bub_idx
Indexes of first & last bubble variable eqns.
real(kind(0d0)) z_centroid
type(bounds_info) z_boundary
These variables combine the centroid and length parameters associated with a particular patch to yiel...
integer gamma_idx
Index of specific heat ratio func. eqn.
type(bounds_info) x_boundary
real(kind(0d0)) eta
In the case that smoothing of patch boundaries is enabled and the boundary between two adjacent patch...
integer, dimension(:,:,:), allocatable logic_grid
procedure(s_assign_patch_xxxxx_primitive_variables), pointer s_assign_patch_primitive_variables
integer e_idx
Index of total energy equation.
subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l)
This subroutine assigns the species primitive variables of the patch designated by the patch_id...
real(kind(0d0)) pref
Reference parameters for Tait EOS.
type(int_bounds_info) mom_idx
Indexes of first & last momentum eqns.
real(kind(0d0)) y_centroid
subroutine s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf)
Converts the conservative variables to the primitive ones.
subroutine s_cylinder(patch_id)
The cylindrical patch is a 3D geometry that may be used, for example, in setting up a cylindrical sol...
This module provides a platform that is analagous to constructive solid geometry techniques and in th...
subroutine s_sweep_plane(patch_id)
The swept plane patch is a 3D geometry that may be used, for example, in creating a solid boundary...
real(kind(0d0)), dimension(:), allocatable pb0
subroutine s_ellipsoid(patch_id)
The ellipsoidal patch is a 3D geometry. The geometry of the patch is well-defined when its centroid a...
subroutine s_3d_analytical(patch_id)
This patch assigns the primitive variables as analytical functions such that the code can be verified...
integer num_patches
Number of patches composing initial condition.
type(ic_patch_parameters), dimension(num_patches_max) patch_icpp
Database of the initial condition patch parameters (icpp) for each of the patches employed in the con...
type(bounds_info) y_boundary
subroutine s_finalize_initial_condition_module()
Deallocation procedures for the module.
real(kind(0d0)), dimension(:), allocatable r0
subroutine s_initialize_initial_condition_module()
Computation of parameters, allocation procedures, and/or any other tasks needed to properly setup the...
subroutine s_generate_initial_condition()
This subroutine peruses the patches and depending on the type of geometry associated with a particula...
type(scalar_field), dimension(:), allocatable q_prim_vf
primitive variables
type(int_bounds_info) internalenergies_idx
Indexes of first & last internal energy eqns.
type(scalar_field) alf_sum
subroutine s_perturb_surrounding_flow()
real(kind(0d0)) sph_phi
Variables to be used to hold cell locations in Cartesian coordinates if 3D simulation is using cylind...
subroutine s_line_segment(patch_id)
The line segment patch is a 1D geometry that may be used, for example, in creating a Riemann problem...
real(kind(0d0)), parameter pi
This module consists of subroutines used in the conversion of the conservative variables into the pri...
logical adv_alphan
Advection of the last volume fraction.
Derived type adding beginning (beg) and end bounds info as attributes.
subroutine s_assign_patch_mixture_primitive_variables(patch_id, j, k, l)
This subroutine assigns the mixture primitive variables of the patch designated by the patch_id...
subroutine s_spherical_harmonic(patch_id)
This patch generates the shape of the spherical harmonics as a perturbation to a perfect sphere...
subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l)
This subroutine assigns the species primitive variables. This follows s_assign_patch_species_primitiv...
subroutine s_3dvarcircle(patch_id)
subroutine s_sweep_line(patch_id)
The swept line patch is a 2D geometry that may be used, for example, in creating a solid boundary...
integer perturb_flow_fluid
Fluid to be perturbed with perturb_flow flag.