[com1] Add additional input options for rheological models

avaframe/com1DFA/DFAfunctionsCython.pyx

@@ -90,17 +90,22 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double sizeSediment = cfg.getfloat('sizeSediment')
   cdef double cvMaxSediment = cfg.getfloat('cvMaxSediment')
   cdef double cvSediment = cfg.getfloat('cvSediment')
+  cdef double etaObrienAndJulienCfg = cfg.getfloat('etaObrienAndJulien')
   cdef double alpha1EtaObrienAndJulien = cfg.getfloat('alpha1EtaObrienAndJulien')
   cdef double beta1EtaObrienAndJulien = cfg.getfloat('beta1EtaObrienAndJulien')
+  cdef double tauyObrienAndJulienCfg = cfg.getfloat('tauyObrienAndJulien')
   cdef double alpha2TauyObrienAndJulien = cfg.getfloat('alpha2TauyObrienAndJulien')
   cdef double beta2TauyObrienAndJulien = cfg.getfloat('beta2TauyObrienAndJulien')
   cdef double alphaObrienAndJulien = cfg.getfloat('alphaObrienAndJulien')
+  cdef double tauyHerschelAndBulkleyCfg = cfg.getfloat('tauyHerschelAndBulkley')
   cdef double alpha2TauyHerschelAndBulkley = cfg.getfloat('alpha2TauyHerschelAndBulkley')
   cdef double beta2TauyHerschelAndBulkley = cfg.getfloat('beta2TauyHerschelAndBulkley')
   cdef double kHerschelAndBulkley = cfg.getfloat('kHerschelAndBulkley')
   cdef double nHerschelAndBulkley = cfg.getfloat('nHerschelAndBulkley')
+  cdef double etaBinghamCfg = cfg.getfloat('etaBingham')
   cdef double alpha1EtaBingham = cfg.getfloat('alpha1EtaBingham')
   cdef double beta1EtaBingham = cfg.getfloat('beta1EtaBingham')
+  cdef double tauyBinghamCfg = cfg.getfloat('tauyBingham')
   cdef double alpha2TauyBingham = cfg.getfloat('alpha2TauyBingham')
   cdef double beta2TauyBingham = cfg.getfloat('beta2TauyBingham')
   cdef double curvAccInFriction = cfg.getfloat('curvAccInFriction')
@@ -168,6 +173,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
   cdef double nx, ny, nz, nxEnd, nyEnd, nzEnd, nxAvg, nyAvg, nzAvg
   cdef double gravAccNorm, accNormCurv, effAccNorm, gravAccTangX, gravAccTangY, gravAccTangZ, forceBotTang, sigmaB, tau
   cdef double muVoellmyRaster, xsiVoellmyRaster
+  #TODO: Influence on variable initiation?
   cdef double shearRate, etaObrienAndJulien, tauyObrienAndJulien, lmObrienAndJulien
   cdef double lambdaBagnold, cObrienAndJulien, tauyHerschelAndBulkley, etaBingham, tauyBingham                                                                          
   # variables for interpolation
@@ -326,9 +332,15 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
           elif frictType == 10:
             ## O`Brien and Julien
             # viscosity
-            etaObrienAndJulien = alpha1EtaObrienAndJulien * math.exp(beta1EtaObrienAndJulien * cvSediment)
+            if etaObrienAndJulienCfg == 0:
+              etaObrienAndJulien = alpha1EtaObrienAndJulien * math.exp(beta1EtaObrienAndJulien * cvSediment)
+            else:
+              etaObrienAndJulien = etaObrienAndJulienCfg
             # yield shear stress
-            tauyObrienAndJulien = alpha2TauyObrienAndJulien * math.exp(beta2TauyObrienAndJulien * cvSediment)
+            if tauyObrienAndJulienCfg == 0:
+              tauyObrienAndJulien = alpha2TauyObrienAndJulien * math.exp(beta2TauyObrienAndJulien * cvSediment)
+            else:
+              tauyObrienAndJulien = tauyObrienAndJulienCfg
             # Prandtl mixing length
             lmObrienAndJulien = 0.4 * h
             # grain concentration
@@ -340,15 +352,24 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType
           elif frictType == 11:
             ## Herschel and Bulkley
             # yield shear stress
-            tauyHerschelAndBulkley = alpha2TauyHerschelAndBulkley * math.exp(beta2TauyHerschelAndBulkley * cvSediment)
+            if tauyHerschelAndBulkleyCfg == 0:
+              tauyHerschelAndBulkley = alpha2TauyHerschelAndBulkley * math.exp(beta2TauyHerschelAndBulkley * cvSediment)
+            else:
+              tauyHerschelAndBulkley = tauyHerschelAndBulkleyCfg
             # shear stress
             tau = tauyHerschelAndBulkley + kHerschelAndBulkley * math.pow(shearRate, nHerschelAndBulkley)
           elif frictType == 12:
             ## Bingham
             # viscosity
-            etaBingham = alpha1EtaBingham * math.exp(beta1EtaBingham * cvSediment)
+            if etaBinghamCfg == 0:
+              etaBingham = alpha1EtaBingham * math.exp(beta1EtaBingham * cvSediment)
+            else:
+              etaBingham = etaBinghamCfg
             # yield shear stress
-            tauyBingham = alpha2TauyBingham * math.exp(beta2TauyBingham * cvSediment)
+            if tauyBinghamCfg == 0:
+              tauyBingham = alpha2TauyBingham * math.exp(beta2TauyBingham * cvSediment)
+            else:
+              tauyBingham = tauyBinghamCfg
             # shear stress
             tau = tauyBingham + etaBingham * shearRate
           else:
@@ -754,7 +775,7 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0):
   cdef double[:] mStoppedArray = np.empty(0, dtype=np.float64)
   cdef double[:] idStoppedArray = np.empty(0, dtype=np.float64)
   cdef double[:] uMagStoppedArray = np.empty(0, dtype=np.float64)
-  cdef long[:] indRemoveParticle
+  cdef long long[:] indRemoveParticle
   # read fields
   cdef double[:] forceX = force['forceX']
   cdef double[:] forceY = force['forceY']

avaframe/com1DFA/checkCfg.py

@@ -99,7 +99,7 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
     Returns
     --------
     cfg: dict
-        upated configuration settings
+        updated configuration settings
     """
 
     frictParameters = [
@@ -130,17 +130,22 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
         "tau0coulombminshear",
         "muvoellmyminshear",
         "xsivoellmyminshear",
+        "etaObrienAndJulien",
         "alpha1EtaObrienAndJulien",
         "beta1EtaObrienAndJulien",
+        "tauyObrienAndJulien",
         "alpha2TauyObrienAndJulien",
         "beta2TauyObrienAndJulien",
         "alphaObrienAndJulien",
+        "tauyHerschelAndBulkley",
         "alpha2TauyHerschelAndBulkley",
         "beta2TauyHerschelAndBulkley",
         "kHerschelAndBulkley",
         "nHerschelAndBulkley",
+        "etaBingham",
         "alpha1EtaBingham",
         "beta1EtaBingham",
+        "tauyBingham",
         "alpha2TauyBingham",
         "beta2TauyBingham",
     ]
@@ -182,20 +187,51 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
                 for frictP in frictParameters
             ]
         elif frictModel.lower() == "obrienandjulien":
-            frictParameterIn = [
-                frictModel.lower() in frictP.lower()
-                for frictP in frictParameters
-            ]
+            if cfg["GENERAL"]["etaObrienAndJulien"] == "0":
+                frictParameters.remove("etaObrienAndJulien")
+            else:
+                frictParameters.remove("alpha1EtaObrienAndJulien")
+                frictParameters.remove("beta1EtaObrienAndJulien")
+                cfg["GENERAL"]["alpha1EtaObrienAndJulien"] = "0"
+                cfg["GENERAL"]["beta1EtaObrienAndJulien"] = "0"
+
+            if cfg["GENERAL"]["tauyObrienAndJulien"] == "0":
+                frictParameters.remove("tauyObrienAndJulien")
+            else:
+                frictParameters.remove("alpha2TauyObrienAndJulien")
+                frictParameters.remove("beta2TauyObrienAndJulien")
+                cfg["GENERAL"]["alpha2TauyObrienAndJulien"] = "0"
+                cfg["GENERAL"]["beta2TauyObrienAndJulien"] = "0"
+
+            frictParameterIn = [frictModel.lower() in frictP.lower() for frictP in frictParameters]
         elif frictModel.lower() == "herschelandbulkley":
-            frictParameterIn = [
-                frictModel.lower() in frictP.lower()
-                for frictP in frictParameters
-            ]
+            if cfg["GENERAL"]["tauyHerschelAndBulkley"] == "0":
+                frictParameters.remove("tauyHerschelAndBulkley")
+            else:
+                frictParameters.remove("alpha2TauyHerschelAndBulkley")
+                frictParameters.remove("beta2TauyHerschelAndBulkley")
+                cfg["GENERAL"]["alpha2TauyHerschelAndBulkley"] = "0"
+                cfg["GENERAL"]["beta2TauyHerschelAndBulkley"] = "0"
+
+            frictParameterIn = [frictModel.lower() in frictP.lower() for frictP in frictParameters]
         elif frictModel.lower() == "bingham":
-            frictParameterIn = [
-                frictModel.lower() in frictP.lower()
-                for frictP in frictParameters
-            ]
+            if cfg["GENERAL"]["etaBingham"] == "0":
+                frictParameters.remove("etaBingham")
+            else:
+                frictParameters.remove("alpha1EtaBingham")
+                frictParameters.remove("beta1EtaBingham")
+                cfg["GENERAL"]["alpha1EtaBingham"] = "0"
+                cfg["GENERAL"]["beta1EtaBingham"] = "0"
+
+            if cfg["GENERAL"]["tauyBingham"] == "0":
+                frictParameters.remove("tauyBingham")
+            else:
+                frictParameters.remove("alpha2TauyBingham")
+                frictParameters.remove("beta2TauyBingham")
+                cfg["GENERAL"]["alpha2TauyBingham"] = "0"
+                cfg["GENERAL"]["beta2TauyBingham"] = "0"
+
+            frictParameterIn = [frictModel.lower() in frictP.lower() for frictP in frictParameters]
         else:
             frictParameterIn = [frictModel.lower() in frictP for frictP in frictParameters]
 
@@ -248,8 +284,8 @@ def checkCfgFrictionModel(cfg, inputSimFiles, relVolume=""):
     # O´Brien and Julien friction model
     # fetch parameters and check if cvSediment < cvMaxSediment
     if frictModel.lower() == "obrienandjulien":
-        cvSediment = cfg['GENERAL']['cvSediment']
-        cvMaxSediment = cfg['GENERAL']['cvMaxSediment']
+        cvSediment = cfg["GENERAL"]["cvSediment"]
+        cvMaxSediment = cfg["GENERAL"]["cvMaxSediment"]
         if cvSediment >= cvMaxSediment:
             message = "cvSediment must be smaller than cvMaxSediment"
             log.error(message)

avaframe/com1DFA/com1DFA.py

@@ -922,8 +922,10 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
             "alpha": cfgGen["alphaObrienAndJulien"],
             "Cmax": cfgGen["cvMaxSediment"],
             "Cv": cfgGen["cvSediment"],
+            "eta": cfgGen["etaObrienAndJulien"],
             "alpha1": cfgGen["alpha1EtaObrienAndJulien"],
             "beta1": cfgGen["beta1EtaObrienAndJulien"],
+            "tauy": cfgGen["tauyObrienAndJulien"],
             "alpha2": cfgGen["alpha2TauyObrienAndJulien"],
             "beta2": cfgGen["beta2TauyObrienAndJulien"],
         }
@@ -933,15 +935,19 @@ def createReportDict(avaDir, logName, relName, inputSimLines, cfg, reportAreaInf
             "model": "Herschel and Bulkley",
             "n": cfgGen["nHerschelAndBulkley"],
             "K": cfgGen["kHerschelAndBulkley"],
+            "tauy": cfgGen["tauyHerschelAndBulkley"],
             "alpha2": cfgGen["alpha2TauyHerschelAndBulkley"],
             "beta2": cfgGen["beta2TauyHerschelAndBulkley"],
         }
     elif cfgGen["frictModel"].lower() == "bingham":
         reportST["Friction model"] = {
             "type": "columns",
             "model": "Bingham",
+            "Cv": cfgGen["cvSediment"],
+            "eta": cfgGen["etaBingham"],
             "alpha1": cfgGen["alpha1EtaBingham"],
             "beta1": cfgGen["beta1EtaBingham"],
+            "tauy": cfgGen["tauyBingham"],
             "alpha2": cfgGen["alpha2TauyBingham"],
             "beta2": cfgGen["beta2TauyBingham"],
         }
@@ -3634,7 +3640,12 @@ def fetchRelVolume(releaseFile, cfg, pathToDem, secondaryReleaseFile, radius=0.0
 
     if cfg["GENERAL"].getboolean("timeDependentRelease"):
         relVolume = initializeRelVol(
-            cfg, demVol, releaseFile, radius, releaseType="timeDepRel", timeDepRelFile=timeDepRelFile
+            cfg,
+            demVol,
+            releaseFile,
+            radius,
+            releaseType="timeDepRel",
+            timeDepRelFile=timeDepRelFile,
         )
     else:
         # compute volume of release area

avaframe/com1DFA/com1DFACfg.ini

@@ -389,39 +389,58 @@ cvSediment = 0.5
 #++++++++++++++Parameter Rheological models
 #++++++++++++O´Brien and Julien
 # viscosity parameter
+## TODO adjusting, default setting
+etaObrienAndJulien = 20
 ## TODO adjusting, default value O´Brien and Julien (1985)
+# set etaObrienAndJulien to 0
 alpha1EtaObrienAndJulien = 0.650
 ## TODO adjusting, default value O´Brien and Julien (1985)
+# set etaObrienAndJulien to 0
 beta1EtaObrienAndJulien = 16.81
 # yield shear stress parameter
+## TODO adjusting, default setting
+tauyObrienAndJulien = 200
 ## TODO adjusting, default value O´Brien and Julien (1985)
+# set tauyObrienAndJulien to 0
 alpha2TauyObrienAndJulien = 0.00886
- ## TODO adjusting, default value O´Brien and Julien (1985)
+## TODO adjusting, default value O´Brien and Julien (1985)
+# set tauyObrienAndJulien to 0
 beta2TauyObrienAndJulien = 13.11
 # parameter, default value Takahashi (1978)
 alphaObrienAndJulien = 0.01
 #++++++++++++Herschel and Bulkley
 # yield shear stress parameter
- ## TODO adjusting, default value O´Brien and Julien (1985)
+tauyHerschelAndBulkley = 200
+## TODO adjusting, default value O´Brien and Julien (1985)
+# set tauyHerschelAndBulkley to 0
 alpha2TauyHerschelAndBulkley = 0.00886
- ## TODO adjusting, default value O´Brien and Julien (1985)
+# set tauyHerschelAndBulkley to 0
+## TODO adjusting, default value O´Brien and Julien (1985)
 beta2TauyHerschelAndBulkley = 13.11
 # consistency factor K
 ## TODO adjusting
-kHerschelAndBulkley = 50
+kHerschelAndBulkley = 20
 # flow exponent
 ## TODO adjusting
 nHerschelAndBulkley = 1.5
 #++++++++++++Bingham
 # viscosity parameter
+## TODO adjusting, default setting
+etaBingham = 20
 ## TODO adjusting, default value O´Brien and Julien (1985)
+# set etaBingham to 0
 alpha1EtaBingham = 0.650
- ## TODO adjusting, default value O´Brien and Julien (1985)
+## TODO adjusting, default value O´Brien and Julien (1985)
+# set etaBingham to 0
 beta1EtaBingham = 16.81
 # yield shear stress parameter
- ## TODO adjusting, default value O´Brien and Julien (1985)
+## TODO adjusting, default setting
+tauyBingham = 200
+## TODO adjusting, default value O´Brien and Julien (1985)
+# set tauyBingham to 0
 alpha2TauyBingham = 0.00886
- ## TODO adjusting, default value O´Brien and Julien (1985)
+## TODO adjusting, default value O´Brien and Julien (1985)
+# set tauyBingham to 0
 beta2TauyBingham = 13.11
 
 #++++++++++++ Resistance model