Merge branch 'antidunes-force-density' into 'master'

Use force density in antidunes showcase

See merge request walberla/walberla!613

apps/showcases/Antidunes/Antidunes.prm

@@ -11,10 +11,10 @@ BlockForestParameters
 DomainParameters
 {
    domainSize         <3200, 60, 160>;
-   wavePeriods        1;      // never set to 0 -> division by zero, even if you initialize a flat particle bed
-   liquidHeightFactor 2.9655; // h_0 / d (water height / avg particle diameter) -> from experiment [E1 = 2.9655, E4 = 3.5862]
-   floorHeightFactor  4.1393;  // from domain bottom to sine's average
-   initialAmplitude   0;     // defined from minimum peak to maximum peak as by Pascal et al. (2021)
+   wavePeriods        1;         // never set to 0 -> division by zero, even if you initialize a flat particle bed
+   liquidHeightFactor 2.862;    // h_0 / d (water height / avg particle diameter) -> from experiment [E1=2.862, E2=3.1724, E3=3.27586, E4=3.5862]
+   floorHeightFactor  4.1393;    // from domain bottom to sine's average
+   initialAmplitude   0;         // defined from minimum peak to maximum peak as by Pascal et al. (2021)
 }
 
 PIDParameters
@@ -32,9 +32,9 @@ PhysicsParameters
 {
    enableWetting     false;
    timesteps         2000000;
-   Re                3100; // [E1=3100, E4=4800]
-   Fr                1.31; // [E1=1.31, E4=1.45]
-   We                15.6188;   // [E1=15.6188 , E4=30.2493]
+   Re                3100;       // [E1=3100, E2=3772, E3=4180, E4=4800]
+   Fr                1.31;       // [E1=1.31, E2=1.38, E3=1.44, E4=1.45]
+   We                15.6188;    // [E1=15.6188, E2=21.48, E3=25.54, E4=30.2493]
 }
 
 ParticleParameters
@@ -55,9 +55,8 @@ ModelParameters
 {
    pdfReconstructionModel        OnlyMissing;
    pdfRefillingModel             EquilibriumRefilling;
-   excessMassDistributionModel   EvenlyAllInterface;
+   excessMassDistributionModel   EvenlyNewInterfaceFallbackLiquid;
    curvatureModel                FiniteDifferenceMethod;
-   enableForceWeighting          false;
    useSimpleMassExchange         false;
    cellConversionThreshold       1e-2;
    cellConversionForceThreshold  1e-1;

apps/showcases/Antidunes/AntidunesLatticeModelGeneration.py

 import sympy as sp
-
+import pystencils as ps
 from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, create_lb_collision_rule
 from lbmpy.enums import ForceModel, Method, Stencil
 from lbmpy.stencils import LBStencil
@@ -7,64 +7,30 @@ from lbmpy.stencils import LBStencil
 from pystencils_walberla import CodeGeneration
 from lbmpy_walberla import generate_lattice_model
 
-# general parameters
-generatedMethod = "Cumulants"  # "TRT", "SRT"
-stencilStr = "D3Q27"
-stencil = LBStencil(Stencil.D3Q19 if stencilStr == "D3Q19" else Stencil.D3Q27)
-force = sp.symbols("force_:3")
-layout = "fzyx"
+with CodeGeneration() as ctx:
+    # general parameters
+    layout = 'fzyx'
+    data_type = "float64" if ctx.double_accuracy else "float32"
+
+    stencil = LBStencil(Stencil.D3Q27)
+    omega = sp.Symbol('omega')
+    force_field = ps.fields(f"force(3): {data_type}[3D]", layout=layout)
 
-if generatedMethod == "Cumulants":
-    omega = sp.Symbol("omega")
-    # method definition
-    lbm_config = LBMConfig(
-        stencil=stencil,
-        method=Method.CUMULANT,
-        relaxation_rate=omega,
-        compressible=True,
-        force=force,
-        zero_centered=False,
-        streaming_pattern="pull",
-        galilean_correction=True if stencil == LBStencil(Stencil.D3Q27) else False,
-    )  # free surface implementation only works with pull pattern
-elif generatedMethod == "TRT":
-    omega_e = sp.Symbol("omega_e")
-    omega_o = sp.Symbol("omega_o")
-    # method definition
-    lbm_config = LBMConfig(
-        stencil=stencil,
-        method=Method.TRT,
-        smagorinsky=False,
-        relaxation_rates=[omega_e, omega_o],
-        compressible=True,
-        force=force,
-        force_model=ForceModel.GUO,
-        zero_centered=False,
-        streaming_pattern="pull",
-    )  # free surface implementation only works with pull pattern
-elif generatedMethod == "SRT":
-    omega = sp.Symbol("omega")
     # method definition
-    lbm_config = LBMConfig(
-        stencil=stencil,
-        method=Method.SRT,
-        smagorinsky=True,
-        relaxation_rate=omega,
-        compressible=True,
-        force=force,
-        force_model=ForceModel.GUO,
-        zero_centered=False,
-        streaming_pattern="pull",
-    )  # free surface implementation only works with pull pattern
+    lbm_config = LBMConfig(stencil=stencil,
+                           method=Method.CUMULANT,
+                           relaxation_rate=omega,
+                           compressible=True,
+                           force=force_field,
+                           zero_centered=False,
+                           streaming_pattern='pull',    # free surface implementation only works with pull pattern
+                           galilean_correction=True)
 
-# optimizations to be used by the code generator
-lbm_opt = LBMOptimisation(cse_global=True, field_layout=layout)
+    # optimizations to be used by the code generator
+    lbm_opt = LBMOptimisation(cse_global=True,
+                              field_layout=layout)
 
-collision_rule = create_lb_collision_rule(
-    lbm_config=lbm_config, lbm_optimisation=lbm_opt
-)
+    collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
+                                              lbm_optimisation=lbm_opt)
 
-with CodeGeneration() as ctx:
-    generate_lattice_model(
-        ctx, "AntidunesLatticeModel", collision_rule, field_layout=layout
-    )
+    generate_lattice_model(ctx, "AntidunesLatticeModel", collision_rule, field_layout=layout)