MFC/src/common/include/1dHardcodedIC.fpp at be79d1b8109832b7661d16287b18ffca3e8d0c73 · MFlowCode/MFC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#:def Hardcoded1DVariables()
    ! Place any declaration of intermediate variables here
    real(wp) :: x_mid_diffu, width_sq, profile_shape, temp, molar_mass_inv, y1, y2, y3, y4
#:enddef

#:def Hardcoded1D()
    select case (patch_icpp(patch_id)%hcid)
    case (150)  ! 1D Smooth Alfven Case for MHD
        ! velocity
        q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x%cc(i))
        q_prim_vf(eqn_idx%mom%beg + 2)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x%cc(i))

        ! magnetic field
        q_prim_vf(eqn_idx%B%end - 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x%cc(i))
        q_prim_vf(eqn_idx%B%end)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x%cc(i))
    case (170)  ! 1D profile from external data (e.g. Cantera, SDtoolbox)
        ! This hardcoded case can be used to start a simulation with initial conditions given from a known 1D profile (e.g. Cantera,
        ! SDtoolbox)
        @: HardcodedReadValues()
    case (180)  ! Shu-Osher problem
        ! This is patch is hard-coded for test suite optimization used in the 1D_shuoser cases: "patch_icpp(2)%alpha_rho(1)": "1 +
        ! 0.2*sin(5*x)"
        if (patch_id == 2) then
            q_prim_vf(eqn_idx%cont%beg + 0)%sf(i, 0, 0) = 1 + 0.2*sin(5*x%cc(i))
        end if
    case (181)  ! Titarev-Torro problem
        ! This is patch is hard-coded for test suite optimization used in the 1D_titarevtorro cases: "patch_icpp(2)%alpha_rho(1)":
        ! "1 + 0.1*sin(20*x*pi)"
        q_prim_vf(eqn_idx%cont%beg + 0)%sf(i, 0, 0) = 1 + 0.1*sin(20*x%cc(i)*pi)
    case (182)  ! Multi-component diffusion
        ! This patch is a hard-coded for test suite optimization (multiple component diffusion)
        x_mid_diffu = 0.05_wp/2.0_wp
        width_sq = (2.5_wp*10.0_wp**(-3.0_wp))**2
        profile_shape = 1.0_wp - 0.5_wp*exp(-(x%cc(i) - x_mid_diffu)**2/width_sq)
        q_prim_vf(eqn_idx%mom%beg)%sf(i, 0, 0) = 0.0_wp
        q_prim_vf(eqn_idx%E)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5
        q_prim_vf(eqn_idx%adv%beg)%sf(i, 0, 0) = 1.0_wp

        y1 = (0.195_wp - 0.142_wp)*profile_shape + 0.142_wp
        y2 = (0.0_wp - 0.1_wp)*profile_shape + 0.1_wp
        y3 = (0.214_wp - 0.0_wp)*profile_shape + 0.0_wp
        y4 = (0.591_wp - 0.758_wp)*profile_shape + 0.758_wp

        q_prim_vf(eqn_idx%species%beg)%sf(i, 0, 0) = y1
        q_prim_vf(eqn_idx%species%beg + 1)%sf(i, 0, 0) = y2
        q_prim_vf(eqn_idx%species%beg + 2)%sf(i, 0, 0) = y3
        q_prim_vf(eqn_idx%species%beg + 3)%sf(i, 0, 0) = y4

        temp = (320.0_wp - 1350.0_wp)*profile_shape + 1350.0_wp

        molar_mass_inv = y1/31.998_wp + y2/18.01508_wp + y3/16.04256_wp + y4/28.0134_wp

        q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)

    case(191)  ! 1D Dual Isothermal case

        q_prim_vf(eqn_idx%E)%sf(i, 0, 0) = 101325.0_wp
        q_prim_vf(eqn_idx%mom%beg)%sf(i, 0, 0) = 0.0_wp
        q_prim_vf(eqn_idx%species%beg)%sf(i, 0, 0) = 1.0_wp

        if (x%cc(i) <= 0.025_wp) then
            temp = 700.0_wp + ((1000.0_wp - 700.0_wp)/0.025_wp)*x%cc(i)
        else
            temp = 1200.0_wp + ((900.0_wp - 1000.0_wp)/0.025_wp)*(x%cc(i) - 0.025_wp)
        end if

        molar_mass_inv = 1.0_wp/2.01588_wp
        q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 101325.0_wp/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
    case default
        call s_int_to_str(patch_id, iStr)
        call s_mpi_abort("Invalid hcid specified for patch " // trim(iStr))
    end select
#:enddef