⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 vb_simple.bas

📁 VB版本_SIMPLE算法和11各个算例
💻 BAS
📖 第 1 页 / 共 2 页
字号:
Attribute VB_Name = "VB_SIMPLE"
Public LStop As Boolean
Public Title(13) As String
Public FLOW#, DIFF#, ACOF#
Public LSolve(10) As Boolean
Public LPrint(13) As Boolean
Public LBlk(10) As Boolean
Public f(22, 22, 10) As Double
Public P(22, 22) As Double
Public RHO(22, 22) As Double
Public GAM(22, 22) As Double
Public CON(22, 22) As Double
Public AIP(22, 22) As Double
Public AIM(22, 22) As Double
Public AJP(22, 22) As Double
Public AJM(22, 22) As Double
Public AP(22, 22) As Double
Public X(22) As Double
Public XU(22) As Double
Public XDIF(22) As Double
Public XCV(22) As Double
Public XCVS(22) As Double
Public Y(22) As Double
Public YV(22) As Double
Public YDIF(22) As Double
Public YCV(22) As Double
Public YCVS(22) As Double
Public YCVR(22) As Double
Public YCVRS(22) As Double
Public ARX(22) As Double
Public ARXJ(22) As Double
Public R(22) As Double
Public RMN(22) As Double
Public SX(22) As Double
Public SXMN(22) As Double
Public XCVI(22) As Double
Public XCVIP(22) As Double
Public DU(22, 22) As Double
Public DV(22, 22) As Double
Public FV(22) As Double
Public FVP(22) As Double
Public FX(22) As Double
Public FXM(22) As Double
Public FY(22) As Double
Public FYM(22) As Double
Public PT(22) As Double
Public QT(22) As Double
Public NF%, NFMAX%, NP%, NRHO%
Public NGAM%, L1%, L2%, L3%
Public M1%, M2%, M3%, IST%
Public JST%, ITER%, LAST%
Public RELAX(13) As Double
Public TIME#, DT#, XL#, YL#, IPREF%, JPREF%
Public MODE%, NTIMES(10) As Integer
Public RHOCON As Double
Public SMAX As Double, SSUM As Double
Public ARXJP(22) As Double
Public IEND%, JEND%

Sub Main()
  If FileExists("计算结果.txt") Then
        Kill "计算结果.txt"
  End If
  Open "计算结果.txt" For Append As #1
  Call Setup0
  Call GRID
  Call setup1
  Call START
  Do
    Call DENSE
    Call BOUND
    Call OUTPUT
    If LStop Then Exit Do
    Call Setup2
  Loop
  Close #1
End Sub
'-----------辅助函数测试文件存在------------------
Private Function FileExists(ByVal f$) As Boolean
    On Error Resume Next
    SetAttr f$, vbNormal
    If Err Then
        FileExists = False
    Else
        FileExists = True
    End If
    Err.Clear
End Function
Sub DIFLOW()
  ACOF = DIFF
  If FLOW = 0# Then
    Exit Sub
 End If
    temp = DIFF - Abs(FLOW) * 0.1
    ACOF = 0#
    If temp <= 0 Then
       Exit Sub
    End If
    temp = temp / DIFF
    ACOF = DIFF * temp ^ 5
End Sub
Function Max(a1 As Double, a2 As Double) As Double
    If a1 > a2 Then
      Max = a1
    Else
      Max = a2
    End If
End Function
Function Min(a1 As Integer, a2 As Integer) As Integer
   If a1 < a2 Then
     Min = a1
   Else
     Min = a2
   End If
End Function
Sub Solve()
      ISTF = IST - 1
      JSTF = JST - 1
      IT1 = L2 + IST
      IT2 = L3 + IST
      JT1 = M2 + JST
      JT2 = M3 + JST
      For nt = 1 To NTIMES(NF)
           N = NF
           If LBlk(NF) Then
               PT(ISTF) = 0#
               QT(ISTF) = 0#
               For i = IST To L2
                     BL = 0#: BLP = 0#: BLM = 0#: BLC = 0#
                     For j = JST To M2
                          BL = BL + AP(i, j)
                          If j <> M2 Then
                              BL = BL - AJP(i, j)
                          End If
                          If j <> JST Then
                              BL = BL - AJM(i, j)
                          End If
                          BLP = BLP + AIP(i, j): BLM = BLM + AIM(i, j)
                          BLC = BLC + CON(i, j) + AIP(i, j) * f(i + 1, j, N) + AIM(i, j) * f(i - 1, j, N) _
                                + AJP(i, j) * f(i, j + 1, N) + AJM(i, j) * f(i, j - 1, N) - AP(i, j) * f(i, j, N)
                     Next j
                     DENOM = BL - PT(i - 1) * BLM: DENO = 1E+15
                     If Abs(DENOM / BL) < 0.0000000001 Then
                          DENOM = 1E+20 * DENO
                     End If
                     PT(i) = BLP / DENOM: QT(i) = (BLC + BLM * QT(i - 1)) / DENOM
               Next i
               BL = 0#
               For ii = IST To L2
                     i = IT1 - ii: BL = BL * PT(i) + QT(i)
                     For j = JST To M2
                          f(i, j, N) = f(i, j, N) + BL
                     Next j
               Next ii
               PT(JSTF) = 0#: QT(JSTF) = 0#
               For j = JST To M2
                     BL = 0#: BLP = 0#: BLM = 0#: BLC = 0#
                     For i = IST To L2
                         BL = BL + AP(i, j)
                         If i <> L2 Then
                             BL = BL - AIP(i, j)
                         End If
                         If i <> IST Then
                             BL = BL - AIM(i, j)
                         End If
                         BLP = BLP + AJP(i, j): BLM = BLM + AJM(i, j)
                         BLC = BLC + CON(i, j) + AIP(i, j) * f(i + 1, j, N) + AIM(i, j) * f(i - 1, j, N) _
                              + AJP(i, j) * f(i, j + 1, N) + AJM(i, j) * f(i, j - 1, N) - AP(i, j) * f(i, j, N)
                     Next i
                     DENOM = BL - PT(j - 1) * BLM
                     If Abs(DENOM / BL) < 0.0000000001 Then
                          DENOM = 1E+20 * DENO
                     End If
                     PT(j) = BLP / DENOM: QT(j) = (BLC + BLM * QT(j - 1)) / DENOM
               Next j
               BL = 0#
               For jj = JST To M2
                   j = JT1 - jj
                   BL = BL * PT(j) + QT(j)
                   For i = IST To L2
                       f(i, j, N) = f(i, j, N) + BL
                   Next i
               Next jj
           End If
           For j = JST To M2
               PT(ISTF) = 0#: QT(ISTF) = f(ISTF, j, N)
               For i = IST To L2
                    DENOM = AP(i, j) - PT(i - 1) * AIM(i, j): PT(i) = AIP(i, j) / DENOM
                    temp = CON(i, j) + AJP(i, j) * f(i, j + 1, N) + AJM(i, j) * f(i, j - 1, N)
                    QT(i) = (temp + AIM(i, j) * QT(i - 1)) / DENOM
               Next i
               For ii = IST To L2
                   i = IT1 - ii: f(i, j, N) = f(i + 1, j, N) * PT(i) + QT(i)
               Next ii
           Next j
           For jj = JST To M3
                j = JT2 - jj: PT(ISTF) = 0#: QT(ISTF) = f(ISTF, j, N)
                For i = IST To L2
                     DENOM = AP(i, j) - PT(i - 1) * AIM(i, j): PT(i) = AIP(i, j) / DENOM
                     temp = CON(i, j) + AJP(i, j) * f(i, j + 1, N) + AJM(i, j) * f(i, j - 1, N)
                     QT(i) = (temp + AIM(i, j) * QT(i - 1)) / DENOM
                Next i
                For ii = IST To L2
                      i = IT1 - ii
                      f(i, j, N) = f(i + 1, j, N) * PT(i) + QT(i)
                Next ii
           Next jj
           For i = IST To L2
                PT(JSTF) = 0#: QT(JSTF) = f(i, JSTF, N)
                For j = JST To M2
                     DENOM = AP(i, j) - PT(j - 1) * AJM(i, j): PT(j) = AJP(i, j) / DENOM
                     temp = CON(i, j) + AIP(i, j) * f(i + 1, j, N) + AIM(i, j) * f(i - 1, j, N)
                     QT(j) = (temp + AJM(i, j) * QT(j - 1)) / DENOM
                Next j
                For jj = JST To M2
                    j = JT1 - jj: f(i, j, N) = f(i, j + 1, N) * PT(j) + QT(j)
                Next jj
           Next i
           For ii = IST To L3
                i = IT2 - ii: PT(JSTF) = 0#: QT(JSTF) = f(i, JSTF, N)
                For j = JST To M2
                    DENOM = AP(i, j) - PT(j - 1) * AJM(i, j): PT(j) = AJP(i, j) / DENOM
                    temp = CON(i, j) + AIP(i, j) * f(i + 1, j, N) + AIM(i, j) * f(i - 1, j, N)
                    QT(j) = (temp + AJM(i, j) * QT(j - 1)) / DENOM
                Next j
                For jj = JST To M2
                       j = JT1 - jj: f(i, j, N) = f(i, j + 1, N) * PT(j) + QT(j)
                Next jj
           Next ii
      Next nt
      For j = 2 To M2
          For i = 2 To L2
              CON(i, j) = 0#: AP(i, j) = 0#
          Next i
      Next j
End Sub

'Have some problems need to be solved
Sub Setup0()
      'Public U(22, 22) As Double, V(22, 22) As Double, PC(22, 22) As Double
      NFMAX = 10: NP = 11: NRHO = 12: NGAM = 13: LStop = False
      For i = 1 To 10
          LSolve(i) = False: LBlk(i) = True: NTIMES(i) = 1
      Next i
      For i = 1 To 13
          LPrint(i) = False: RELAX(i) = 1#
      Next i
      MODE = 1: LAST = 5: TIME = 0#
      ITER = 0: DT = 10000000000#: IPREF = 1: JPREF = 1: RHOCON = 1
End Sub
Sub setup1()
      Dim str As String
      L2 = L1 - 1: L3 = L2 - 1
      M2 = M1 - 1: M3 = M2 - 1: X(1) = XU(2)
      For i = 2 To L2
          X(i) = 0.5 * (XU(i + 1) + XU(i))
      Next i
      X(L1) = XU(L1): Y(1) = YV(2)
      For j = 2 To M2
          Y(j) = 0.5 * (YV(j + 1) + YV(j))
      Next j
      Y(M1) = YV(M1)
      For i = 2 To L1
          XDIF(i) = X(i) - X(i - 1)
      Next i
      For i = 2 To L2
          XCV(i) = XU(i + 1) - XU(i)
      Next
      For i = 3 To L2
          XCVS(i) = XDIF(i)
      Next
      XCVS(3) = XCVS(3) + XDIF(2): XCVS(L2) = XCVS(L2) + XDIF(L1)
      For i = 3 To L3
          XCVI(i) = 0.5 * XCV(i)
          XCVIP(i) = XCVI(i)
      Next
      XCVIP(2) = XCV(2)
      XCVI(L2) = XCV(L2)
      For j = 2 To M1
          YDIF(j) = Y(j) - Y(j - 1)
      Next
      For j = 2 To M2
          YCV(j) = YV(j + 1) - YV(j)
      Next
      For j = 3 To M2
          YCVS(j) = YDIF(j)
      Next
      YCVS(3) = YCVS(3) + YDIF(2)
      YCVS(M2) = YCVS(M2) + YDIF(M1)
      If MODE = 1 Then
         For j = 1 To M1
             RMN(j) = 1#
             R(j) = 1#
         Next
      Else
         For j = 2 To M1
             R(j) = R(j - 1) + YDIF(j)
         Next
         RMN(2) = R(1)
         For j = 3 To M2
             RMN(j) = RMN(j - 1) + YCV(j - 1)
         Next
         RMN(M1) = R(M1)
      End If
      
      For j = 1 To M1
          SX(j) = 1#: SXMN(j) = 1#
          If MODE = 3 Then
             SX(j) = R(j)
             If j <> 1 Then
                SXMN(j) = RMN(j)
             End If
          End If
      Next j
      For j = 2 To M2
           YCVR(j) = R(j) * YCV(j): ARX(j) = YCVR(j)
           If MODE = 3 Then
                 ARX(j) = YCV(j)
           End If
      Next j
      For j = 4 To M3
          YCVRS(j) = 0.5 * (R(j) + R(j - 1)) * YDIF(j)
      Next
      YCVRS(3) = 0.5 * (R(3) + R(1)) * YCVS(3)
      YCVRS(M2) = 0.5 * (R(M1) + R(M3)) * YCVS(M2)
      
      If MODE = 2 Then
          For j = 3 To M3
               ARXJ(j) = 0.25 * (1# + RMN(j) / R(j)) * ARX(j)
               ARXJP(j) = ARX(j) - ARXJ(j)
          Next j
      Else
          For j = 3 To M3
              ARXJ(j) = 0.5 * ARX(j): ARXJP(j) = ARXJ(j)
          Next j
      End If
      
      ARXJP(2) = ARX(2): ARXJ(M2) = ARX(M2)
      For j = 3 To M3
          FV(j) = ARXJP(j) / ARX(j)
          FVP(j) = 1# - FV(j)
      Next j
      For i = 3 To L2
          FX(i) = 0.5 * XCV(i - 1) / XDIF(i)
          FXM(i) = 1# - FX(i)
      Next i
      FX(2) = 0#: FXM(2) = 1#: FX(L1) = 1#: FXM(L1) = 0#
      For j = 3 To M2
          FY(j) = 0.5 * YCV(j - 1) / YDIF(j)
          FYM(j) = 1# - FY(j)
      Next
      FY(2) = 0#: FYM(2) = 1#: FY(M1) = 1#: FYM(M1) = 0#
'---CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE----
      For j = 1 To M1
          For i = 1 To L1
              f(i, j, 3) = 0#: f(i, j, 1) = 0#:
              f(i, j, 2) = 0#: CON(i, j) = 0#: AP(i, j) = 0#
              RHO(i, j) = RHOCON: P(i, j) = 0#
          Next i
      Next j
       Print #1,
       Print #1,
       str = Space(15)
      If MODE = 1 Then
         Print #1, str + "Computeation in cartisian coordinates"
      ElseIf MODE = 2 Then
         Print #1, str + "Computation for axisymmetrical situation"
      ElseIf MODE = 3 Then
         Print #1, str + "Computation polar coordinates "
      End If
      Print #1, str + String(40, "*")
End Sub
Sub Setup2()
'---COEFFICIENTS FOR THE U EQUATION----
      NF = 1
      If LSolve(NF) Then
       IST = 3
       JST = 2
       Call GAMSOR
       REL = 1# - RELAX(NF)
       For i = 3 To L2
          FL = XCVI(i) * f(i, 2, 2) * RHO(i, 1)
          FLM = XCVIP(i - 1) * f(i - 1, 2, 2) * RHO(i - 1, 1)
          FLOW = R(1) * (FL + FLM)
          DIFF = R(1) * (XCVI(i) * GAM(i, 1) + XCVIP(i - 1) * GAM(i - 1, 1)) / YDIF(2)
          Call DIFLOW
          AJM(i, 2) = ACOF + Max(0#, FLOW)
       Next i
       For j = 2 To M2
          FLOW = ARX(j) * f(2, j, 1) * RHO(1, j)
          DIFF = ARX(j) * GAM(1, j) / (XCV(2) * SX(j))
          Call DIFLOW
          AIM(3, j) = ACOF + Max(0#, FLOW)
          For i = 3 To L2
              If i <> L2 Then
                 FL = f(i, j, 1) * (FX(i) * RHO(i, j) + FXM(i) * RHO(i - 1, j))
                 FLP = f(i + 1, j, 1) * (FX(i + 1) * RHO(i + 1, j) + FXM(i + 1) * RHO(i, j))
                 FLOW = ARX(j) * 0.5 * (FL + FLP)
                 DIFF = ARX(j) * GAM(i, j) / (XCV(i) * SX(j))
              Else
                 FLOW = ARX(j) * f(L1, j, 1) * RHO(L1, j)

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -