Convergence speeded up

TODO

@@ -1,2 +0,0 @@
-* fvSolution: cellLimited
-* collapseDict: collapseEdges

TODO.md

@@ -0,0 +1,8 @@
+* fvSolution: cellLimited
+* collapseDict: collapseEdges
+
+## 1.03.21
+* boundary type (wall or symetryPlane)
+* restruct for ways
+* build alpha = 0.01 .. 0.13
+* ! symetryPlane -> cyclicAMI

notes.md

@@ -0,0 +1,28 @@
+## relaxationFactors
+
+| p     | Iterations    | U     |
+| ---   | ---           | ---   |
+| 0.05  | 830           | 0.2   |
+| 0.1   | 678           | 0.2   |
+| 0.2   | 505           | 0.2   |
+| 0.3   | 305           | 0.3   |
+| 0.3   | 236           | 0.4   |
+| 0.3   | 181           | 0.5   |
+
+-->
+
+## SIMPLE.residualControl
+
+| p     | U     | Iterations    |
+| ---   | ---   | ---           |
+| 1e-4  | 1e-4  | 338           |
+| 1e-5  | 1e-5  | 499           |
+
+##
+
+U = \frac{\Delta p L}{mu} \frac{9}{256} (\sqrt 2 - \frac{1}{1 - \alpha})^2
+L = 2 R_0
+R_0 = 1
+scale = 1e-5
+mu = 1e-3
+\Delta p = 1

src/genmesh.py

@@ -10,8 +10,12 @@ if not os.path.exists(build):
 
 ###
 
-alpha = [0.1] #, 0.15, 0.2]
+alpha = []
 
+for n in range(0.01, 0.13 + 0.01, 0.01):
+    alpha.append(n)
+
+# Structures
 simpleCubic = os.path.join(src, "simple-cubic/main.py")
 # Body-centered cubic
 #bcCubic = os.path.join(path, "bc-cubic/main.py")
@@ -22,24 +26,28 @@ simpleCubic = os.path.join(src, "simple-cubic/main.py")
 
 processes = []
 structure = ["simple-cubic"] #, "bc-cubic", "fc-cubic"]
+direction = ["1"]
 
 def salome(src_path, build_path, arg):
     subprocess.run(["salome", "start", "-t", src_path, "args:{},{}".format(build_path, arg)])
 
 for s in structure:
     s_path = os.path.join(build, s)
+    
+    for d in direction:
+        d_path = os.path.join(s_path, d)
 
-    for c in alpha:
+        for c in alpha:
-        src_path = os.path.join(src, "%s/main.py" % s)
+            src_path = os.path.join(src, "%s/main.py" % s)
-        build_path = os.path.join(s_path, str(c))
+            build_path = os.path.join(d_path, str(c))
-        
+            
-        if not os.path.exists(build_path):
+            if not os.path.exists(build_path):
-            os.makedirs(build_path)
+                os.makedirs(build_path)
 
-        print("starting process")
+            print("starting process")
-        p = multiprocessing.Process(target = salome, args = (src_path, build_path, c))
+            p = multiprocessing.Process(target = salome, args = (src_path, build_path, c))
-        processes.append(p)
+            processes.append(p)
-        p.start()
+            p.start()
 
-    for process in processes:
+        for process in processes:
-        process.join()
+            process.join()

src/simple-cubic/0/U

@@ -25,7 +25,7 @@ boundaryField
     inlet
     {
         type            fixedValue;
-        value           uniform (0 0 -0.001);
+        value           uniform (0 0 -6e-5);
     }
     outlet
     {

src/simple-cubic/geometry.py

@@ -41,7 +41,14 @@ def create(alpha):
     # Main geometry
     Pore = geompy.MakeCutList(cut3, [sphere4], True)
     geompy.addToStudy(Pore, "PORE")
-   
+
+    #
+    #   D /\ B  A(1, 1, 0)  \vec{n}(1, 1, 0)
+    #    /  \   B(3, 3, 0)  \vec{n}(1, 1, 0)
+    #    \  /   C(3, 1, 0)  \vec{n}(-1, 1, 0)
+    #   A \/ C  D(1, 3, 0)  \vec{n}(-1, 1, 0)
+    #
+
     # Prepare faces
     vtx.append(geompy.MakeVertex(2, 2, 2))
     vtx.append(geompy.MakeVertexWithRef(vtx[3], 0, 0, 1))
@@ -51,6 +58,10 @@ def create(alpha):
     plane = geompy.MakeTranslation(plane, 0, 0, -2)
     common2 = geompy.MakeCommonList([Pore, plane], True)
     
+    # TODO: make 4 planes A B C D
+    Pore = geompy.MakeFillet(Pore, 0.1, geompy.ShapeType["EDGE"], [24, 27, 35])
+
+
     # Main groups (inlet, outlet, wall)
     inlet = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"])
     gip = geompy.GetInPlace(Pore, common1, True)

src/simple-cubic/system/fvSolution

@@ -75,8 +75,8 @@ SIMPLE
 
     residualControl
     {
-        p               1e-3;
+        p               1e-5;
-        U               1e-3;
+        U               1e-5;
     }
 }
 
@@ -84,11 +84,11 @@ relaxationFactors
 {
     fields
     {
-        p               0.05;
+        p               0.3;
     }
     equations
     {
-        U               0.2;
+        U               0.5;
     }
 }
 

src/tools/run.sh

@@ -10,8 +10,8 @@ path=.
 rm -rf postProcessing processor* *.log logs *.obj constant/polyMesh
 
 # Export and transform mesh
-ideasUnvFoam -case $path mesh.unv
+ideasUnvToFoam -case $path mesh.unv
-transformPoints -scale '(0.001 0.001 0.001)'
+transformPoints -scale '(1e-5 1e-5 1e-5)'
 #polyDualMesh 70 -overwrite
 checkMesh -case $path -allGeometry -allTopology > checkMesh.log
 
@@ -22,7 +22,7 @@ foamDictionary -case $path constant/polyMesh/boundary -entry entry0.wall.type -s
 decomposePar -case $path
 
 # Initial approximation
-potentialFoam -case $path -parallel
+mpirun -np 4 --oversubscribe potentialFoam -case $path -parallel
 
 # Change boundary type for simpleFoam
 for n in {0..3}; do
@@ -31,6 +31,7 @@ for n in {0..3}; do
     foamDictionary "processor${n}/0/U" -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
 done
 
+sleep 2
 # Main calculation
-mpirun -np 4 --oversubscribe simpleFoam -parallel -case $path > $path/simpleFoam.log
+mpirun -np 4 --oversubscribe simpleFoam -parallel -case $path | tee -a $path/simpleFoam.log
 

worksheet/cube_010320.py

@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+
+###
+### This file is generated automatically by SALOME v9.6.0 with dump python functionality
+###
+
+import sys
+import salome
+
+salome.salome_init()
+import salome_notebook
+notebook = salome_notebook.NoteBook()
+sys.path.insert(0, r'/home/nafaryus/projects/anisotrope-cube/worksheet')
+
+###
+### GEOM component
+###
+
+import GEOM
+from salome.geom import geomBuilder
+import math
+import SALOMEDS
+
+
+geompy = geomBuilder.New()
+
+geomObj_1 = geompy.MakeVectorDXDYDZ(1, 0, 0)
+geomObj_2 = geompy.MakeVectorDXDYDZ(0, 1, 0)
+geomObj_3 = geompy.MakeVectorDXDYDZ(0, 0, 1)
+geomObj_4 = geompy.MakeBoxDXDYDZ(2.82842712474619, 2.82842712474619, 2)
+geomObj_5 = geompy.MakeRotation(geomObj_4, geomObj_3, 45*math.pi/180.0)
+geomObj_6 = geompy.MakeTranslation(geomObj_5, 2, 0, 0)
+geomObj_7 = geompy.MakeVertex(2, 0, 0)
+geomObj_8 = geompy.MakeVertex(2, 2, 0)
+geomObj_9 = geompy.MakeVertex(2, 2, 2)
+geomObj_10 = geompy.MakeLineTwoPnt(geomObj_8, geomObj_9)
+geomObj_11 = geompy.MakeSpherePntR(geomObj_7, 1.111111111111111)
+geomObj_12 = geompy.MakeMultiTranslation1D(geomObj_11, geomObj_2, 2, 3)
+geomObj_13 = geompy.MakeCutList(geomObj_6, [geomObj_12], True)
+geomObj_14 = geompy.MakeTranslation(geomObj_12, 0, 0, 2)
+geomObj_15 = geompy.MakeCutList(geomObj_13, [geomObj_14], True)
+geomObj_16 = geompy.MakeRotation(geomObj_12, geomObj_10, 90*math.pi/180.0)
+geomObj_17 = geompy.MakeCutList(geomObj_15, [geomObj_16], True)
+geomObj_18 = geompy.MakeTranslation(geomObj_16, 0, 0, 2)
+PORE = geompy.MakeCutList(geomObj_17, [geomObj_18], True)
+geomObj_19 = geompy.MakeVertex(2, 2, 2)
+geomObj_20 = geompy.MakeVertexWithRef(geomObj_19, 0, 0, 1)
+geomObj_21 = geompy.MakeVector(geomObj_19, geomObj_20)
+geomObj_22 = geompy.MakePlane(geomObj_19, geomObj_21, 5)
+geomObj_23 = geompy.MakeCommonList([PORE, geomObj_22], True)
+geomObj_24 = geompy.MakeTranslation(geomObj_22, 0, 0, -2)
+geomObj_25 = geompy.MakeCommonList([PORE, geomObj_24], True)
+
+inlet = geompy.CreateGroup(PORE, geompy.ShapeType["FACE"])
+geomObj_26 = geompy.GetInPlace(PORE, geomObj_23, True)
+[geomObj_27,geomObj_28,geomObj_29,geomObj_30] = geompy.SubShapeAll(geomObj_26, geompy.ShapeType["FACE"])
+geompy.UnionList(inlet, [geomObj_27, geomObj_28, geomObj_29, geomObj_30])
+outlet = geompy.CreateGroup(PORE, geompy.ShapeType["FACE"])
+geomObj_31 = geompy.GetInPlace(PORE, geomObj_25, True)
+[geomObj_32,geomObj_33,geomObj_34,geomObj_35] = geompy.SubShapeAll(geomObj_31, geompy.ShapeType["FACE"])
+geompy.UnionList(outlet, [geomObj_32, geomObj_33, geomObj_34, geomObj_35])
+
+geomObj_36 = geompy.CreateGroup(PORE, geompy.ShapeType["FACE"])
+geompy.UnionIDs(geomObj_36, [3, 17, 28, 33, 42, 51, 56, 63, 72, 77, 90, 95, 102, 107, 117, 120, 129, 136, 142, 149, 156, 159, 162, 165])
+wall = geompy.CutListOfGroups([geomObj_36], [inlet, outlet])
+
+Fillet_1 = geompy.MakeFillet(PORE, 0.1, geompy.ShapeType["EDGE"], [24, 27, 35, 41, 48, 58, 82, 85, 89, 109, 114, 122, 128, 135, 140, 141, 148, 155])
+
+inlet_1 = geompy.CreateGroup(Fillet_1, geompy.ShapeType["FACE"])
+geompy.UnionIDs(inlet_1, [138, 37, 96, 269])
+outlet_1 = geompy.CreateGroup(Fillet_1, geompy.ShapeType["FACE"])
+geompy.UnionIDs(outlet_1, [64, 172, 126, 262])
+wall_1 = geompy.CutListOfGroups([geomObj_37], [inlet_1, outlet_1])
+geompy.addToStudy( PORE, 'PORE' )
+geompy.addToStudyInFather( PORE, inlet, 'inlet' )
+geompy.addToStudyInFather( PORE, outlet, 'outlet' )
+geompy.addToStudyInFather( PORE, wall, 'wall' )
+geompy.addToStudy( Fillet_1, 'Fillet_1' )
+geompy.addToStudyInFather( Fillet_1, inlet_1, 'inlet' )
+geompy.addToStudyInFather( Fillet_1, outlet_1, 'outlet' )
+geompy.addToStudyInFather( Fillet_1, wall_1, 'wall' )
+
+###
+### SMESH component
+###
+
+import  SMESH, SALOMEDS
+from salome.smesh import smeshBuilder
+
+smesh = smeshBuilder.New()
+#smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed or in some particular situations:
+                                 # multiples meshes built in parallel, complex and numerous mesh edition (performance)
+
+PORE_1 = smesh.Mesh(PORE)
+NETGEN_2D3D_1 = PORE_1.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
+NETGEN_Parameters = NETGEN_2D3D_1.Parameters()
+NETGEN_Parameters.SetSecondOrder( 0 )
+NETGEN_Parameters.SetOptimize( 1 )
+NETGEN_Parameters.SetChordalError( -1 )
+NETGEN_Parameters.SetChordalErrorEnabled( 0 )
+NETGEN_Parameters.SetUseSurfaceCurvature( 1 )
+NETGEN_Parameters.SetFuseEdges( 1 )
+NETGEN_Parameters.SetCheckChartBoundary( 0 )
+NETGEN_Parameters.SetMinSize( 0.01 )
+NETGEN_Parameters.SetMaxSize( 0.1 )
+NETGEN_Parameters.SetFineness( 4 )
+NETGEN_Parameters.SetQuadAllowed( 0 )
+ViscousLayers_0_025 = NETGEN_2D3D_1.ViscousLayers(0.025,2,1.1,[],1,smeshBuilder.NODE_OFFSET)
+inlet_2 = PORE_1.GroupOnGeom(inlet,'inlet',SMESH.FACE)
+outlet_2 = PORE_1.GroupOnGeom(outlet,'outlet',SMESH.FACE)
+wall_2 = PORE_1.GroupOnGeom(wall,'wall',SMESH.FACE)
+isDone = PORE_1.Compute()
+try:
+  PORE_1.ExportUNV( r'asd/mesh.unv' )
+  pass
+except:
+  print('ExportUNV() failed. Invalid file name?')
+NETGEN_Parameters.SetFineness( 3 )
+Mesh_1 = smesh.Mesh(Fillet_1)
+status = Mesh_1.AddHypothesis(NETGEN_Parameters)
+status = Mesh_1.AddHypothesis(ViscousLayers_0_025)
+NETGEN_2D3D_1_1 = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
+inlet_3 = Mesh_1.GroupOnGeom(inlet_1,'inlet',SMESH.FACE)
+outlet_3 = Mesh_1.GroupOnGeom(outlet_1,'outlet',SMESH.FACE)
+wall_3 = Mesh_1.GroupOnGeom(wall_1,'wall',SMESH.FACE)
+ViscousLayers_0_025.SetTotalThickness( 0.05 )
+ViscousLayers_0_025.SetNumberLayers( 2 )
+ViscousLayers_0_025.SetStretchFactor( 1.1 )
+ViscousLayers_0_025.SetMethod( smeshBuilder.SURF_OFFSET_SMOOTH )
+ViscousLayers_0_025.SetFaces( [], 1 )
+isDone = Mesh_1.Compute()
+[ inlet_3, outlet_3, wall_3 ] = Mesh_1.GetGroups()
+try:
+  Mesh_1.ExportUNV( r'/home/nafaryus/projects/anisotrope-cube/build/simple-cubic/0.1/mesh.unv' )
+  pass
+except:
+  print('ExportUNV() failed. Invalid file name?')
+
+
+## Set names of Mesh objects
+smesh.SetName(NETGEN_2D3D_1.GetAlgorithm(), 'NETGEN_2D3D_1')
+smesh.SetName(ViscousLayers_0_025, 'ViscousLayers=0.025,2,1.1,[],1')
+smesh.SetName(NETGEN_Parameters, 'NETGEN_Parameters')
+smesh.SetName(inlet_2, 'inlet')
+smesh.SetName(outlet_2, 'outlet')
+smesh.SetName(wall_2, 'wall')
+smesh.SetName(PORE_1.GetMesh(), 'PORE')
+smesh.SetName(Mesh_1.GetMesh(), 'Mesh_1')
+smesh.SetName(wall_3, 'wall')
+smesh.SetName(outlet_3, 'outlet')
+smesh.SetName(inlet_3, 'inlet')
+
+
+if salome.sg.hasDesktop():
+  salome.sg.updateObjBrowser()