Mod: track new model variables

Mod: return mesh elements
Mod: update netgen conda channel
Fix: some mistakes

.vscode/launch.json

@@ -0,0 +1,15 @@
+{
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "name": "anisotropy compute",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/anisotropy/core/cli.py",
+            "args": ["compute", "-v"],
+            "console": "integratedTerminal",
+            "cwd": "/tmp/anisotropy",
+            "preLaunchTask": "prepare"
+        }
+    ]
+}

.vscode/settings.json

@@ -0,0 +1,12 @@
+{
+    "python.pythonPath": "${env:HOME}/.local/share/mambaforge/envs/anisotropy/bin/python",
+    "python.testing.pytestArgs": [
+        "tests"
+    ],
+    "python.testing.unittestEnabled": false,
+    "python.testing.pytestEnabled": true,
+    "python.linting.enabled": true,
+    "python.linting.pylintEnabled": false,
+    "python.linting.flake8Enabled": true,
+    "python.linting.flake8Args": ["--ignore=E402,E251,E501", "--verbose"]
+}

.vscode/tasks.json

@@ -0,0 +1,9 @@
+{
+    "version": "2.0.0",
+    "tasks": [{
+        "label": "prepare",
+        "command": "mkdir",
+        "args": ["-p", "/tmp/anisotropy"],
+        "type": "shell",
+    }]
+}

anisotropy/core/cli.py

@@ -4,7 +4,7 @@
 
 import click
 import ast
-import os, sys, shutil
+import os
 import logging
 
 

anisotropy/core/config.py

@@ -72,22 +72,22 @@ class Config(object):
         #   Expand structures for each direction and each alpha
         for structure in self.content["structures"]:
             # ISSUE: precision error 0.06999999999999999
+            structure = deepcopy(structure)
             alphaA = round(arange(
                 structure["alpha"][0], 
                 structure["alpha"][1] + structure["alphaStep"], 
                 structure["alphaStep"]
             ), 9)
             directionA = array(structure["directions"], dtype = float) 
+            structure.pop("alpha")
+            structure.pop("directions")
 
             for direction in directionA:
                 for alpha in alphaA:
                     self.cases.append({
-                        "label": structure["label"],
                         "alpha": alpha,
-                        "alphaStep": structure["alphaStep"],
                         "direction": direction,
-                        "r0": structure["r0"],
+                        **structure
-                        "filletsEnabled": structure["filletsEnabled"]
                     })
 
     def chooseParams(self, idn: int):
@@ -122,7 +122,15 @@ class DefaultConfig(Config):
                 "label": label,
                 "alpha": alpha,
                 "alphaStep": 0.01,
-                "directions": [[1, 0, 0], [0, 0, 1], [1, 1, 1]],
+                "directions": [[1., 0., 0.], [0., 0., 1.], [1., 1., 1.]],
                 "r0": 1,
-                "filletsEnabled": True
+                "filletsEnabled": True,
+                "pressureInlet": 1,
+                "pressureOutlet": 0,
+                "pressureInternal": 0,
+                "velocityInlet": [0., 0., 0.],
+                "velocityOutlet": None,
+                "velocityInternal": [0., 0., 0.],
+                "density": 1000,
+                "viscosity": 1e-3
             })

anisotropy/core/postProcess.py

@@ -16,6 +16,7 @@ class PostProcess(object):
         self.path = path.abspath(dirpath)
 
     def flowRate(self, patch: str):
+        # TODO: fix wrong log path
         func = "patchFlowRate(patch={})".format(patch)
         filepath = path.join(self.path, "postProcessing", func, "0", "surfaceFieldValue.dat")
         postProcess(func, cwd = self.path)

anisotropy/core/runner.py

@@ -31,6 +31,7 @@ class UltimateRunner(object):
 
         # Database preparation
         self.database = Database(path = self.config["database"])
+        self.exec_id = None
 
         if exec_id:
             if self.database.getExecution(exec_id):
@@ -64,7 +65,7 @@ class UltimateRunner(object):
             flow = self.database.getFlowOnephase(execution = self.exec_id, **self.config.params)
 
             if not flow:
-                flow = T.FlowOnephase.create(mesh_id = mesh)
+                flow = T.FlowOnephase.create(mesh_id = mesh, **self.config.params)
 
     def fill(self):
         self.config.expand()
@@ -202,6 +203,9 @@ class UltimateRunner(object):
 
         if not returncode:
             meshParams.meshStatus = "done"
+            meshParams.edges = len(self.mesh.edges)
+            meshParams.faces = len(self.mesh.faces)
+            meshParams.volumes = len(self.mesh.volumes)
 
         else:
             logger.error(err)
@@ -225,7 +229,16 @@ class UltimateRunner(object):
         ))
         timer = Timer()
 
-        self.flow = OnePhaseFlow(params["direction"], path = self.casepath())
+        flowParams.viscosityKinematic = flowParams.viscosity / flowParams.density
+
+        with self.database:
+            flowParams.save()
+
+            self.flow = OnePhaseFlow(
+                direction = params["direction"], 
+                **flowParams.select().dicts().get(), 
+                path = self.casepath()
+            )
 
         if not self.shape:
             filename = "shape.step"
@@ -300,8 +313,6 @@ class UltimateRunner(object):
         if stage in ["postProcess", "all"]:
             self.computePostProcess()
 
-        #logger.info("Pipeline done")
-
     @staticmethod
     def subrunner(*args, **kwargs):
         runner = UltimateRunner(config = kwargs["config"], exec_id = kwargs["exec_id"], typo = "worker")

anisotropy/database/models.py

@@ -11,7 +11,7 @@ from peewee import (
     TimeField, DateTimeField, Proxy
 )
 from .utils import JSONField
-#from playhouse.sqlite_ext import JSONField
+
 __database_proxy__ = Proxy()
 
 class Execution(Model):
@@ -45,9 +45,7 @@ class Shape(Model):
 
     volumeCell = FloatField(null = True)
     volume = FloatField(null = True)
-    volumeRounded = FloatField(null = True)
     porosity = FloatField(null = True)
-    porosityRounded = FloatField(null = True)
 
     class Meta:
         database = __database_proxy__
@@ -87,10 +85,11 @@ class FlowOnephase(Model):
     pressureInlet = FloatField(null = True)
     pressureOutlet = FloatField(null = True)
     pressureInternal = FloatField(null = True)
-    velocityInlet = FloatField(null = True)
+    velocityInlet = JSONField(null = True)
-    velocityOutlet = FloatField(null = True)
+    velocityOutlet = JSONField(null = True)
-    velocityInternal = FloatField(null = True)
+    velocityInternal = JSONField(null = True)
     viscosity = FloatField(null = True)
+    viscosityKinematic = FloatField(null = True)
     density = FloatField(null = True)
     flowRate = FloatField(null = True)
     permeability = FloatField(null = True)

anisotropy/database/utils.py

@@ -15,6 +15,7 @@ from peewee import (
 import json
 from numpy import ndarray
 
+
 class ListField(TextField):
     field_type = "list"
 
@@ -173,3 +174,4 @@ class JSONField(TextField):
 
     def tree(self):
         return fn.json_tree(self)
+

anisotropy/meshing/mesh.py

@@ -4,7 +4,7 @@
 
 from netgen.occ import OCCGeometry
 from netgen import meshing
-import numpy
+from numpy import array
 import os
 
 class Mesh(object):
@@ -82,9 +82,20 @@ class Mesh(object):
     def doubleExport(self):
         pass
 
-    def volumes(self) -> numpy.array:
+    @property
-        #   TODO: check each polyhedron
+    def volumes(self) -> array:
-        tetras = numpy.array([ [ [ vertex for vertex in mesh[index] ] for index in element.vertices ] for element in self.mesh.Elements3D() ])
+        return array(self.mesh.Elements3D())
-        volumes = numpy.array([ 1 / 6 * linalg.det(numpy.append(tetra.transpose(), numpy.array([[1, 1, 1, 1]]), axis = 0)) for tetra in tetras ])
+
+    @property
+    def faces(self) -> array:
+        return array(self.mesh.Elements2D())
+
+    @property
+    def edges(self) -> array:
+        # NOTE: returns zero elements for imported mesh 
+        return array(self.mesh.Elements1D())
+
+
+# tetras = numpy.array([ [ [ vertex for vertex in mesh[index] ] for index in element.vertices ] for element in self.mesh.Elements3D() ])
+# volumes = numpy.array([ 1 / 6 * linalg.det(numpy.append(tetra.transpose(), numpy.array([[1, 1, 1, 1]]), axis = 0)) for tetra in tetras ])
         
-        return volumes

anisotropy/shaping/__init__.py

@@ -2,7 +2,7 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from .shape import Shape
+from .shape import ShapeError, Shape
 from .periodic import Periodic
 from .simple import Simple
 from .faceCentered import FaceCentered

anisotropy/shaping/bodyCentered.py

@@ -8,6 +8,7 @@ from numpy import pi, sqrt, arccos
 
 from .occExtended import *
 from . import Periodic
+from . import ShapeError
 
 class BodyCentered(Periodic):
     def __init__(
@@ -160,6 +161,13 @@ class BodyCentered(Periodic):
         #   Main shape
         self.shape = self.cell - self.lattice
 
+        if not len(self.shape.solids) == 1:
+            raise ShapeError("Expected single solid shape")
+
+        else:
+            self.shape = self.shape.solids[0]
+
+        #   Boundaries (walls)
         for face in self.shape.faces:
             if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
                 face.name = "wall"

anisotropy/shaping/faceCentered.py

@@ -8,6 +8,7 @@ from numpy import pi, sqrt
 
 from .occExtended import *
 from . import Periodic
+from . import ShapeError
 
 class FaceCentered(Periodic):
     def __init__(
@@ -167,6 +168,13 @@ class FaceCentered(Periodic):
         #   Main shape
         self.shape = self.cell - self.lattice
 
+        if not len(self.shape.solids) == 1:
+            raise ShapeError("Expected single solid shape")
+
+        else:
+            self.shape = self.shape.solids[0]
+
+        #   Boundaries (walls)
         for face in self.shape.faces:
             if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
                 face.name = "wall"

anisotropy/shaping/shape.py

@@ -7,6 +7,11 @@ import numpy
 from numpy import linalg
 import os
 
+
+class ShapeError(Exception):
+    pass
+
+
 class Shape(object):
     def __init__(self):
         self.groups = {}

anisotropy/shaping/simple.py

@@ -8,6 +8,7 @@ from numpy import pi, sqrt
 
 from .occExtended import *
 from . import Periodic
+from . import ShapeError
 
 class Simple(Periodic):
     def __init__(
@@ -155,6 +156,13 @@ class Simple(Periodic):
         #   Main shape
         self.shape = self.cell - self.lattice
 
+        if not len(self.shape.solids) == 1:
+            raise ShapeError("Expected single solid shape")
+
+        else:
+            self.shape = self.shape.solids[0]
+
+        #   Boundaries (walls)
         for face in self.shape.faces:
             if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
                 face.name = "wall"

anisotropy/solving/onephase.py

@@ -11,9 +11,32 @@ import logging
 logger = logging.getLogger(__name__)
 
 class OnePhaseFlow(FoamCase):
-    def __init__(self, direction, path: str = None):
+    def __init__(
+        self,
+        direction: list[float] = None,
+        pressureInlet: float = None,
+        pressureOutlet: float = None,
+        pressureInternal: float = None,
+        velocityInlet: float = None,
+        velocityOutlet: float = None,
+        velocityInternal: float = None,
+        density: float = None,
+        viscosityKinematic: float = None,
+        path: str = None,
+        **kwargs
+    ):
         FoamCase.__init__(self, path = path)
 
+        self.direction = direction 
+        self.pressureInlet = pressureInlet 
+        self.pressureOutlet = pressureOutlet 
+        self.pressureInternal = pressureInternal 
+        self.velocityInlet = velocityInlet 
+        self.velocityOutlet = velocityOutlet 
+        self.velocityInternal = velocityInternal 
+        self.density = density 
+        self.viscosityKinematic = viscosityKinematic
+
         controlDict = F.ControlDict()
         controlDict.update(
             startFrom = "latestTime",
@@ -27,7 +50,8 @@ class OnePhaseFlow(FoamCase):
         fvSolution = F.FvSolution()
         fvSolution["solvers"]["U"].update(
             nSweeps = 2,
-            tolerance = 1e-08
+            tolerance = 1e-08,
+            smoother = "GaussSeidel"
         )
         fvSolution["solvers"]["Phi"] = dict(
             solver = "GAMG",
@@ -40,7 +64,7 @@ class OnePhaseFlow(FoamCase):
             relTol = 0.01
         )
         fvSolution["potentialFlow"] = dict(
-            nNonOrthogonalCorrectors = 20,
+            nNonOrthogonalCorrectors = 13,
             PhiRefCell = 0,
             PhiRefPoint = 0,
             PhiRefValue = 0,
@@ -48,7 +72,7 @@ class OnePhaseFlow(FoamCase):
         )
         fvSolution["cache"] = { "grad(U)": None }
         fvSolution["SIMPLE"].update(
-            nNonOrthogonalCorrectors = 10,
+            nNonOrthogonalCorrectors = 6,
             residualControl = dict(
                 p = 1e-05,
                 U = 1e-05
@@ -58,7 +82,7 @@ class OnePhaseFlow(FoamCase):
 
         transportProperties = F.TransportProperties()
         transportProperties.update(
-            nu = 1e-06
+            nu = self.viscosityKinematic
         )
         
         turbulenceProperties = F.TurbulenceProperties()
@@ -72,22 +96,21 @@ class OnePhaseFlow(FoamCase):
         u = F.U()
         u["boundaryField"] = {}
 
-        # ISSUE: add proxy from geometry direction to outlet boundaryField.
         for boundary in boundaries:
             if boundary == "inlet":
                 p["boundaryField"][boundary] = dict(
                     type = "fixedValue",
-                    value = uniform(1e-3)
+                    value = uniform(self.pressureInlet / self.density)
                 )
                 u["boundaryField"][boundary] = dict(
                     type = "fixedValue",
-                    value = uniform(array(direction) * -6e-5) # uniform([0, 0, -6e-5])
+                    value = uniform(array(self.direction) * -6e-5)
                 )
 
             elif boundary == "outlet":
                 p["boundaryField"][boundary] = dict(
                     type = "fixedValue",
-                    value = uniform(0)
+                    value = uniform(self.pressureOutlet / self.density)
                 )
                 u["boundaryField"][boundary] = dict(
                     type = "zeroGradient",
@@ -173,13 +196,13 @@ class OnePhaseFlow(FoamCase):
                 "scale": [1e-5, 1e-5, 1e-5]
             })
             R.renumberMesh()
-            R.potentialFoam()
+            #R.potentialFoam()
 
             self.read()
 
-            self.U["boundaryField"]["outlet"] = dict(
+            self.U["boundaryField"]["inlet"] = dict(
                 type = "pressureInletVelocity",
-                value = uniform([0, 0, 0])
+                value = uniform(self.velocityInlet)
             )
             self.write()
 

environment.yml

@@ -8,4 +8,4 @@ dependencies:
   - poetry
   - sqlite
   - occt
-  - netgen
+  - l-nafaryus::netgen