Refactoring force models

doc/conf.py

@@ -21,6 +21,9 @@ extensions = [
     'sphinx_autodoc_typehints',
 ]
 
+# set mathjax v3 path according to https://www.sphinx-doc.org/en/master/usage/extensions/math.html
+mathjax_path="https://cdn.jsdelivr.net/npm/mathjax@2/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
+
 add_module_names = False
 templates_path = ['_templates']
 source_suffix = '.rst'

doc/sphinx/lbmpy.bib

@@ -6,7 +6,8 @@
   number={4},
   pages={043309},
   year={2015},
-  publisher={APS}
+  publisher={APS},
+  doi={10.1103/PhysRevE.92.043309},
 }
 
 @PHDTHESIS{luo1993lattice,
@@ -16,7 +17,7 @@
    school = {GEORGIA INSTITUTE OF TECHNOLOGY.},
      year = 1993,
    adsurl = {http://adsabs.harvard.edu/abs/1993PhDT.......233L},
-  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System},
 }
 
 @Article{guo2002discrete,
@@ -28,8 +29,24 @@
   number    = {4},
   pages     = {046308},
   publisher = {APS},
+  doi = {10.1103/PhysRevE.65.046308},
 }
 
+@article{HeForce,
+  title = {Discrete Boltzmann equation model for nonideal gases},
+  author = {He, Xiaoyi and Shan, Xiaowen and Doolen, Gary D.},
+  journal = {Physical Review E},
+  volume = {57},
+  issue = {1},
+  pages = {R13--R16},
+  numpages = {0},
+  year = {1998},
+  month = {1},
+  publisher = {APS},
+  doi = {10.1103/PhysRevE.57.R13}
+}
+
+
 @article{buick2000gravity,
   title={Gravity in a lattice Boltzmann model},
   author={Buick, JM and Greated, CA},
@@ -38,7 +55,8 @@
   number={5},
   pages={5307},
   year={2000},
-  publisher={APS}
+  publisher={APS},
+  doi = {10.1103/PhysRevE.61.5307},
 }
 
 @phdthesis{schiller2008thermal,

lbmpy_tests/test_force.py

@@ -9,6 +9,7 @@ from contextlib import ExitStack as does_not_raise
 
 
 force_models = [fm.lower() for fm in dir(lbmpy.forcemodels) if fm[0].isupper()]
+force_models.remove("abstractforcemodel")
 
 def test_force_models_available():
     assert 'guo' in force_models
@@ -82,7 +83,7 @@ def test_total_momentum(method, force_model, omega):
 
 
 @pytest.mark.parametrize("stencil", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
-@pytest.mark.parametrize("force_model", ["simple", "schiller"])
+@pytest.mark.parametrize("force_model", ["simple", "schiller", "Guo"])
 def test_modes(stencil, force_model):
     """check Schiller's force term in mode space"""
     stencil = get_stencil(stencil)
@@ -101,8 +102,11 @@ def test_modes(stencil, force_model):
                               compressible=True,
                               force_model=force_model,
                               force=F)
-    
-    force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
+
+    if force_model == "Guo":
+        force_moments = method.force_model.moment_space_forcing(method)
+    else:
+        force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
     
     # The mass mode should be zero
     assert force_moments[0] == 0
@@ -110,7 +114,7 @@ def test_modes(stencil, force_model):
     # The momentum moments should contain the force
     assert list(force_moments[1:dim+1]) == F
     
-    if force_model == "schiller":
+    if force_model == "schiller" or force_model=="Guo":
         num_stresses = (dim*dim-dim)//2+dim
         lambda_s, lambda_b = -omega_s, -omega_b
 
@@ -134,7 +138,10 @@ def test_modes(stencil, force_model):
             else:
                 assert diff == 0 # difference should be zero
 
-        ff = sp.Matrix(method.force_model(method))
+        if force_model == "Guo":
+            ff = method.moment_matrix.inv() * sp.Matrix(method.force_model.moment_space_forcing(method))
+        else:
+            ff = sp.Matrix(method.force_model(method))
         # Check eq. 4.53a from schiller2008thermal
         assert sp.simplify(sum(ff)) == 0
         # Check eq. 4.53b from schiller2008thermal
@@ -192,3 +199,21 @@ def test_momentum_density_shift(force_model):
     dh.run_kernel(momentum_density_kernel)
     assert np.sum(dh.gather_array(momentum_density.name)[:, :, 0]) == np.prod(domain_size) / 2
     assert np.sum(dh.gather_array(momentum_density.name)[:, :, 1]) == np.prod(domain_size)
+
+
+@pytest.mark.parametrize('fmodel_class', [lbmpy.forcemodels.Simple, lbmpy.forcemodels.Luo, lbmpy.forcemodels.Guo])
+def test_forcing_space_equivalences(fmodel_class):
+    d, q = 3, 27
+    stencil = get_stencil(f"D{d}Q{q}")
+    force = sp.symbols(f"F_:{d}")
+    fmodel = fmodel_class(force)
+    lb_method = create_lb_method(stencil=stencil, method='srt')
+    inv_moment_matrix = lb_method.moment_matrix.inv()
+
+    force_pdfs = sp.Matrix(fmodel(lb_method))
+    force_moments = fmodel.moment_space_forcing(lb_method)
+
+    diff = (force_pdfs - (inv_moment_matrix * force_moments)).expand()
+    for i, d in enumerate(diff):
+        assert d == 0, f"Mismatch between population and moment space forcing " \
+                       f"in force model {fmodel_class}, population f_{i}"
\ No newline at end of file