diff --git a/TODO b/TODO deleted file mode 100644 index 7aa300d..0000000 --- a/TODO +++ /dev/null @@ -1,2 +0,0 @@ -* fvSolution: cellLimited -* collapseDict: collapseEdges diff --git a/TODO.md b/TODO.md new file mode 100644 index 0000000..7a4dceb --- /dev/null +++ b/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 diff --git a/notes.md b/notes.md new file mode 100644 index 0000000..02a1f05 --- /dev/null +++ b/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 diff --git a/src/genmesh.py b/src/genmesh.py index 2314b8d..7af546c 100644 --- a/src/genmesh.py +++ b/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: - src_path = os.path.join(src, "%s/main.py" % s) - build_path = os.path.join(s_path, str(c)) - - if not os.path.exists(build_path): - os.makedirs(build_path) + for c in alpha: + src_path = os.path.join(src, "%s/main.py" % s) + build_path = os.path.join(d_path, str(c)) + + if not os.path.exists(build_path): + os.makedirs(build_path) - print("starting process") - p = multiprocessing.Process(target = salome, args = (src_path, build_path, c)) - processes.append(p) - p.start() + print("starting process") + p = multiprocessing.Process(target = salome, args = (src_path, build_path, c)) + processes.append(p) + p.start() - for process in processes: - process.join() + for process in processes: + process.join() diff --git a/src/simple-cubic/0/U b/src/simple-cubic/0/U index ab060ba..c0e8bf3 100644 --- a/src/simple-cubic/0/U +++ b/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 { diff --git a/src/simple-cubic/geometry.py b/src/simple-cubic/geometry.py index 4e876ad..66241c1 100644 --- a/src/simple-cubic/geometry.py +++ b/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) diff --git a/src/simple-cubic/system/fvSolution b/src/simple-cubic/system/fvSolution index d32674e..09e75c5 100644 --- a/src/simple-cubic/system/fvSolution +++ b/src/simple-cubic/system/fvSolution @@ -75,8 +75,8 @@ SIMPLE residualControl { - p 1e-3; - U 1e-3; + p 1e-5; + U 1e-5; } } @@ -84,11 +84,11 @@ relaxationFactors { fields { - p 0.05; + p 0.3; } equations { - U 0.2; + U 0.5; } } diff --git a/src/tools/run.sh b/src/tools/run.sh index 7a7c206..6170651 100755 --- a/src/tools/run.sh +++ b/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 -transformPoints -scale '(0.001 0.001 0.001)' +ideasUnvToFoam -case $path mesh.unv +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 diff --git a/worksheet/cube_010320.py b/worksheet/cube_010320.py new file mode 100644 index 0000000..cc6d2e1 --- /dev/null +++ b/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()