diff --git a/examples/Hybrid/Mixing1a.JSON b/examples/Hybrid/Mixing1a.JSON new file mode 100644 index 0000000..bf49225 --- /dev/null +++ b/examples/Hybrid/Mixing1a.JSON @@ -0,0 +1,243 @@ +{ +"network": { + "nodes": [ + { + "x": 0.0, + "y": 0.0, + "z": 0.0, + "ground": true + }, + { + "x": 1e-3, + "y": 2e-3, + "z": 0.0 + }, + { + "x": 1e-3, + "y": 1e-3, + "z": 0.0 + }, + { + "x": 1e-3, + "y": 0.0, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 2e-3, + "z": 0.0 + }, + { + "x": 1.75e-3, + "y": 1e-3, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 0.0, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 1.25e-3, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 0.75e-3, + "z": 0.0 + }, + { + "x": 2.25e-3, + "y": 1e-3, + "z": 0.0 + }, + { + "x": 3e-3, + "y": 1e-3, + "z": 0.0, + "ground": true + } + ], + "channels": [ + { + "node1": 0, + "node2": 1, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 0, + "node2": 2, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 0, + "node2": 3, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 1, + "node2": 4, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 2, + "node2": 5, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 3, + "node2": 6, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 4, + "node2": 7, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 6, + "node2": 8, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 9, + "node2": 10, + "width": 1e-4, + "height": 1e-4 + } + ], + "modules": [ + { + "position": [1.75e-3, 0.75e-3], + "size": [5e-4, 5e-4], + "nodes": [5, 7, 8, 9] + } + ] +}, +"simulation": { + "platform": "Mixing", + "type": "Hybrid", + "resistanceModel": "Poiseuille", + "mixingModel": "Instantaneous", + "fluids": [ + { + "name": "Water", + "concentration": 1, + "density": 1000, + "viscosity": 0.001 + } + ], + "species": [ + { + "name": "O2", + "diffusivity": 1e-8, + "saturationConcentration": 1.0, + "molecularSize": 0.0 + } + ], + "mixtures": [ + { + "species": [0], + "concentrations": [1.1] + } + ], + "pumps": [ + { + "channel":0, + "type": "PumpPressure", + "deltaP": 1000 + }, + { + "channel":1, + "type": "PumpPressure", + "deltaP": 1000 + }, + { + "channel":2, + "type": "PumpPressure", + "deltaP": 1000 + } + ], + "fixtures": [ + { + "name": "Setup #1", + "phase": 0, + "mixtureInjections": [ + { + "mixture": 0, + "channel": 4, + "t0": 0.0 + } + ] + } + ], + "activeFixture": 0, + "settings": { + "simulators": [ + { + "Type": "Mixing", + "name": "Mixing1a", + "stlFile": "../examples/STL/cross.stl", + "charPhysLength": 1e-4, + "charPhysVelocity": 1e-1, + "alpha": 0.1, + "resolution": 20, + "epsilon": 1e-1, + "tau": 0.55, + "moduleId": 0, + "Openings": [ + { + "node": 5, + "normal": { + "x": 1.0, + "y": 0.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + }, + { + "node": 7, + "normal": { + "x": 0.0, + "y": -1.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + }, + { + "node": 8, + "normal": { + "x": 0.0, + "y": 1.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + }, + { + "node": 9, + "normal": { + "x": -1.0, + "y": 0.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + } + ] + } + ] + } +} +} \ No newline at end of file diff --git a/examples/Hybrid/OoC-1.JSON b/examples/Hybrid/OoC-1.JSON new file mode 100644 index 0000000..1fb9a56 --- /dev/null +++ b/examples/Hybrid/OoC-1.JSON @@ -0,0 +1,253 @@ +{ +"network": { + "nodes": [ + { + "x": 0.0, + "y": 0.0, + "z": 0.0, + "ground": true + }, + { + "x": 1e-3, + "y": 2e-3, + "z": 0.0 + }, + { + "x": 1e-3, + "y": 1e-3, + "z": 0.0 + }, + { + "x": 1e-3, + "y": 0.0, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 2e-3, + "z": 0.0 + }, + { + "x": 1.75e-3, + "y": 1e-3, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 0.0, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 1.25e-3, + "z": 0.0 + }, + { + "x": 2e-3, + "y": 0.75e-3, + "z": 0.0 + }, + { + "x": 2.25e-3, + "y": 1e-3, + "z": 0.0 + }, + { + "x": 3e-3, + "y": 1e-3, + "z": 0.0, + "ground": true + } + ], + "channels": [ + { + "node1": 0, + "node2": 1, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 0, + "node2": 2, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 0, + "node2": 3, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 1, + "node2": 4, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 2, + "node2": 5, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 3, + "node2": 6, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 4, + "node2": 7, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 6, + "node2": 8, + "width": 1e-4, + "height": 1e-4 + }, + { + "node1": 9, + "node2": 10, + "width": 1e-4, + "height": 1e-4 + } + ], + "modules": [ + { + "position": [1.75e-3, 0.75e-3], + "size": [5e-4, 5e-4], + "nodes": [5, 7, 8, 9] + } + ] +}, +"simulation": { + "platform": "Ooc", + "type": "Hybrid", + "resistanceModel": "Poiseuille", + "mixingModel": "Instantaneous", + "fluids": [ + { + "name": "Water", + "concentration": 1, + "density": 1000, + "viscosity": 0.001 + } + ], + "tissues": [ + { + "name": "Liver", + "species": [0], + "Vmax": [0.1571428], + "kM": [0.022285714] + } + ], + "species": [ + { + "name": "O2", + "diffusivity": 1e-8, + "saturationConcentration": 1.0, + "molecularSize": 0.0 + } + ], + "mixtures": [ + { + "species": [0], + "concentrations": [1.0] + } + ], + "pumps": [ + { + "channel":0, + "type": "PumpPressure", + "deltaP": 1000 + }, + { + "channel":1, + "type": "PumpPressure", + "deltaP": 1000 + }, + { + "channel":2, + "type": "PumpPressure", + "deltaP": 1000 + } + ], + "fixtures": [ + { + "name": "Setup #1", + "phase": 0, + "mixtureInjections": [ + { + "mixture": 0, + "channel": 1, + "t0": 0.0 + } + ] + } + ], + "activeFixture": 0, + "settings": { + "simulators": [ + { + "Type": "Organ", + "name": "Liver-Test", + "stlFile": "../examples/STL/cross.stl", + "tissue": 0, + "organStlFile": "../examples/STL/cube.stl", + "charPhysLength": 1e-4, + "charPhysVelocity": 1e-1, + "alpha": 0.1, + "resolution": 20, + "epsilon": 1e-1, + "tau": 0.55, + "moduleId": 0, + "Openings": [ + { + "node": 5, + "normal": { + "x": 1.0, + "y": 0.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + }, + { + "node": 7, + "normal": { + "x": 0.0, + "y": -1.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + }, + { + "node": 8, + "normal": { + "x": 0.0, + "y": 1.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + }, + { + "node": 9, + "normal": { + "x": -1.0, + "y": 0.0, + "z": 0.0 + }, + "width": 1e-4, + "height": 1e-4 + } + ] + } + ] + } +} +} \ No newline at end of file diff --git a/examples/STL/cube.stl b/examples/STL/cube.stl new file mode 100644 index 0000000..230b92d --- /dev/null +++ b/examples/STL/cube.stl @@ -0,0 +1,198 @@ +solid +facet normal 0.000000e+00 -1.000000e+00 0.000000e+00 + outer loop + vertex 2.250000e-04 2.250000e-04 2.500000e-05 + vertex 2.750000e-04 2.250000e-04 2.500000e-05 + vertex 2.250000e-04 2.250000e-04 7.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 -1.000000e+00 0.000000e+00 + outer loop + vertex 2.750000e-04 2.250000e-04 2.500000e-05 + vertex 2.750000e-04 2.250000e-04 7.500000e-05 + vertex 2.250000e-04 2.250000e-04 7.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 1.000000e+00 0.000000e+00 + outer loop + vertex 2.250000e-04 2.750000e-04 2.500000e-05 + vertex 2.250000e-04 2.750000e-04 7.500000e-05 + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 1.000000e+00 0.000000e+00 + outer loop + vertex 2.250000e-04 2.750000e-04 7.500000e-05 + vertex 2.750000e-04 2.750000e-04 7.500000e-05 + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + endloop +endfacet +facet normal -1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.250000e-04 2.250000e-04 2.500000e-05 + vertex 2.250000e-04 2.250000e-04 7.500000e-05 + vertex 2.250000e-04 2.750000e-04 2.500000e-05 + endloop +endfacet +facet normal -1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.250000e-04 2.250000e-04 7.500000e-05 + vertex 2.250000e-04 2.750000e-04 7.500000e-05 + vertex 2.250000e-04 2.750000e-04 2.500000e-05 + endloop +endfacet +facet normal 1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.750000e-04 2.250000e-04 2.500000e-05 + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + vertex 2.750000e-04 2.250000e-04 7.500000e-05 + endloop +endfacet +facet normal 1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + vertex 2.750000e-04 2.750000e-04 7.500000e-05 + vertex 2.750000e-04 2.250000e-04 7.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 1.000000e+00 + outer loop + vertex 2.250000e-04 2.250000e-04 7.500000e-05 + vertex 2.750000e-04 2.250000e-04 7.500000e-05 + vertex 2.250000e-04 2.750000e-04 7.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 1.000000e+00 + outer loop + vertex 2.750000e-04 2.250000e-04 7.500000e-05 + vertex 2.750000e-04 2.750000e-04 7.500000e-05 + vertex 2.250000e-04 2.750000e-04 7.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 1.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 1.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 -1.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 -1.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal -1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal -1.000000e+00 0.000000e+00 0.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.250000e-04 2.250000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.750000e-04 2.250000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.250000e-04 2.250000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.250000e-04 2.750000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.250000e-04 2.250000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.250000e-04 2.750000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.250000e-04 2.750000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.750000e-04 2.250000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.750000e-04 2.750000e-04 2.500000e-05 + endloop +endfacet +facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 + outer loop + vertex 2.750000e-04 2.250000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + vertex 2.500000e-04 2.500000e-04 2.500000e-05 + endloop +endfacet +endsolid \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1b14291..59a7056 100755 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,6 +10,7 @@ set(HEADER_LIST add_subdirectory(nodalAnalysis) add_subdirectory(architecture) add_subdirectory(simulation) +add_subdirectory(olbProcessors) add_subdirectory(porting) add_subdirectory(result) diff --git a/src/baseSimulator.h b/src/baseSimulator.h index e6698fd..91df00c 100644 --- a/src/baseSimulator.h +++ b/src/baseSimulator.h @@ -7,6 +7,7 @@ #include "simulation/Droplet.h" #include "simulation/Fluid.h" #include "simulation/Injection.h" +#include "simulation/MembraneModels.h" #include "simulation/Mixture.h" #include "simulation/MixtureInjection.h" #include "simulation/MixingModels.h" @@ -14,7 +15,10 @@ #include "simulation/Simulation.h" #include "simulation/simulators/cfdSimulator.h" #include "simulation/simulators/olbContinuous.h" +#include "simulation/simulators/olbMixing.h" +#include "simulation/simulators/olbOoc.h" #include "simulation/Specie.h" +#include "simulation/Tissue.h" #include "simulation/events/BoundaryEvent.h" #include "simulation/events/Event.h" #include "simulation/events/InjectionEvent.h" @@ -32,6 +36,10 @@ #include "architecture/Node.h" #include "architecture/PressurePump.h" +#include "olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.h" +#include "olbProcessors/saturatedFluxPostProcessor2D.h" +#include "olbProcessors/setFunctionalRegularizedHeatFlux.h" + #include "porting/jsonPorter.h" #include "porting/jsonReaders.h" #include "porting/jsonWriters.h" @@ -40,4 +48,5 @@ #ifdef USE_ESSLBM #include "simulation/simulators/essContinuous.h" + #include "simulation/simulators/essMixing.h" #endif \ No newline at end of file diff --git a/src/baseSimulator.hh b/src/baseSimulator.hh index 4aa4c28..b2e5192 100644 --- a/src/baseSimulator.hh +++ b/src/baseSimulator.hh @@ -2,6 +2,7 @@ #include "simulation/Droplet.hh" #include "simulation/Fluid.hh" #include "simulation/Injection.hh" +#include "simulation/MembraneModels.hh" #include "simulation/Mixture.hh" #include "simulation/MixtureInjection.hh" #include "simulation/MixingModels.hh" @@ -9,7 +10,10 @@ #include "simulation/Simulation.hh" #include "simulation/simulators/cfdSimulator.hh" #include "simulation/simulators/olbContinuous.hh" +#include "simulation/simulators/olbMixing.hh" +#include "simulation/simulators/olbOoc.hh" #include "simulation/Specie.hh" +#include "simulation/Tissue.hh" #include "simulation/events/BoundaryEvent.hh" #include "simulation/events/InjectionEvent.hh" #include "simulation/events/MergingEvent.hh" @@ -25,6 +29,10 @@ #include "architecture/Node.hh" #include "architecture/PressurePump.hh" +#include "olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh" +#include "olbProcessors/saturatedFluxPostProcessor2D.hh" +#include "olbProcessors/setFunctionalRegularizedHeatFlux.hh" + #include "porting/jsonPorter.hh" #include "porting/jsonReaders.hh" #include "porting/jsonWriters.hh" @@ -33,4 +41,5 @@ #ifdef USE_ESSLBM #include "simulation/simulators/essContinuous.hh" + #include "simulation/simulators/essMixing.hh" #endif \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 37fffd1..ca4b9f1 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,7 +1,6 @@ //#include //namespace py = pybind11; - #include #include diff --git a/src/nodalAnalysis/NodalAnalysis.hh b/src/nodalAnalysis/NodalAnalysis.hh index 1325beb..5527c92 100644 --- a/src/nodalAnalysis/NodalAnalysis.hh +++ b/src/nodalAnalysis/NodalAnalysis.hh @@ -302,7 +302,7 @@ void NodalAnalysis::initGroundNodes(std::unordered_mapsetGroundNodes(groundNodes); - cfdSimulator->setFlowRates(flowRates_); + cfdSimulator->storeFlowRates(flowRates_); cfdSimulator->setInitialized(true); } } @@ -398,8 +398,8 @@ void NodalAnalysis::writeCfdSimulators(std::unordered_mapsetPressures(pressures_); - cfdSimulator.second->setFlowRates(flowRates_); + cfdSimulator.second->storePressures(pressures_); + cfdSimulator.second->storeFlowRates(flowRates_); } } diff --git a/src/olbProcessors/CMakeLists.txt b/src/olbProcessors/CMakeLists.txt new file mode 100644 index 0000000..f350763 --- /dev/null +++ b/src/olbProcessors/CMakeLists.txt @@ -0,0 +1,15 @@ +set(SOURCE_LIST + navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh + saturatedFluxPostProcessor2D.hh + setFunctionalRegularizedHeatFlux.hh +) + +set(HEADER_LIST + navierStokesAdvectionDiffusionCouplingPostProcessor2D.h + saturatedFluxPostProcessor2D.h + setFunctionalRegularizedHeatFlux.h +) + +target_sources(${TARGET_NAME} PUBLIC ${SOURCE_LIST} ${HEADER_LIST}) +target_include_directories(${TARGET_NAME} PUBLIC ${CMAKE_CURRENT_LIST_DIR}) +target_link_libraries(${TARGET_NAME} PUBLIC lbmLib) \ No newline at end of file diff --git a/src/olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.h b/src/olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.h new file mode 100644 index 0000000..3473426 --- /dev/null +++ b/src/olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.h @@ -0,0 +1,52 @@ +#pragma once + +#include +#include + +namespace olb { + +/** +* Class for the coupling between a Navier-Stokes (NS) lattice and an +* Advection-Diffusion (AD) lattice. +*/ + +//====================================================================== +// ========== AD coupling in Single Direction 2D =======================// +//====================================================================== +template +class NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D : public LocalPostProcessor2D { +private: + typedef DESCRIPTOR L; + int x0, x1, y0, y1; + const std::vector velFactors; + std::vector>*> tPartners; + std::vector*> partners; + +public: + NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D(int x0_, int x1_, int y0_, int y1_, std::vector velFactors, + std::vector* > partners_); + int extent() const override + { + return 0; + } + int extent(int whichDirection) const override + { + return 0; + } + void process(BlockLattice& blockLattice) override; + void processSubDomain(BlockLattice& blockLattice, + int x0_, int x1_, int y0_, int y1_) override; +}; + +template +class NavierStokesAdvectionDiffusionSingleCouplingGenerator2D : public LatticeCouplingGenerator2D { +private: + const std::vector velFactors; + +public: + NavierStokesAdvectionDiffusionSingleCouplingGenerator2D(int x0_, int x1_, int y0_, int y1_, std::vector velFactors_); + PostProcessor2D* generate(std::vector* > partners) const override; + LatticeCouplingGenerator2D* clone() const override; +}; + +} diff --git a/src/olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh b/src/olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh new file mode 100644 index 0000000..4cd98a0 --- /dev/null +++ b/src/olbProcessors/navierStokesAdvectionDiffusionCouplingPostProcessor2D.hh @@ -0,0 +1,78 @@ +#include "navierStokesAdvectionDiffusionCouplingPostProcessor2D.h" + +namespace olb { + +//===================================================================================== +//=========== NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D ============ +//===================================================================================== + +template +NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D:: +NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D(int x0_, int x1_, int y0_, int y1_, std::vector velFactors_, std::vector* > partners_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), velFactors(velFactors_), partners(partners_) +{ + this->getName() = "NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D"; + for (long unsigned int i = 0; i> *>(partners[i])); + } +} + +template +void NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D:: +processSubDomain(BlockLattice& blockLattice, + int x0_, int x1_, int y0_, int y1_) +{ + int newX0, newX1, newY0, newY1; + if ( util::intersect ( + x0, x1, y0, y1, + x0_, x1_, y0_, y1_, + newX0, newX1, newY0, newY1 ) ) { + + for (long unsigned int i = 0; iget(iX,iY); + // Velocity coupling + T u[DESCRIPTOR::d] { }; + cell.computeU(u); + for ( int d = 0; d(u); + } + } + } + } +} + +template +void NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D:: +process(BlockLattice& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1); +} + +/// LatticeCouplingGenerator for advectionDiffusion coupling + +template +NavierStokesAdvectionDiffusionSingleCouplingGenerator2D:: +NavierStokesAdvectionDiffusionSingleCouplingGenerator2D(int x0_, int x1_, int y0_, int y1_, std::vector velFactors_) + : LatticeCouplingGenerator2D(x0_, x1_, y0_, y1_), velFactors(velFactors_) +{ } + +template +PostProcessor2D* NavierStokesAdvectionDiffusionSingleCouplingGenerator2D::generate ( + std::vector* > partners) const +{ + return new NavierStokesAdvectionDiffusionSingleCouplingPostProcessor2D( + this->x0,this->x1,this->y0,this->y1, this->velFactors, partners); +} + +template +LatticeCouplingGenerator2D* NavierStokesAdvectionDiffusionSingleCouplingGenerator2D::clone() const +{ + return new NavierStokesAdvectionDiffusionSingleCouplingGenerator2D(*this); +} + +} // namespace olb diff --git a/src/olbProcessors/saturatedFluxPostProcessor2D.h b/src/olbProcessors/saturatedFluxPostProcessor2D.h new file mode 100644 index 0000000..c4c31ce --- /dev/null +++ b/src/olbProcessors/saturatedFluxPostProcessor2D.h @@ -0,0 +1,47 @@ +#pragma once + +#include +#include + +namespace olb { + +//////// SaturationFluxProcessor2D //////////////////////////// + +template +class SaturatedFluxProcessor2D : public LocalPostProcessor2D { +public: + SaturatedFluxProcessor2D(int x0_, int x1_, int y0_, int y1_, + int discreteNormalX_, int discreteNormalY_, T *Vmax, T *Km); + int extent() const override + { + return 1; + } + int extent(int whichDirection) const override + { + return 1; + } + void process(BlockLattice& blockLattice) override; + void processSubDomain(BlockLattice& blockLattice, + int x0_,int x1_,int y0_,int y1_ ) override; +private: + int x0, x1, y0, y1; + int discreteNormalX, discreteNormalY; + T *Vmax; + T *Km; +}; + +/// Generator class for the FreeEnergyWall PostProcessor handling the wetting boundary condition. +template +class SaturatedFluxProcessorGenerator2D : public PostProcessorGenerator2D { +public: + SaturatedFluxProcessorGenerator2D(int x0_, int x1_, int y0_, int y1_, + int discreteNormalX_, int discreteNormalY_, T *Vmax_, T *Km_); + PostProcessor2D* generate() const override; + PostProcessorGenerator2D* clone() const override; +private: + int discreteNormalX; + int discreteNormalY; + T *Vmax; + T *Km; +}; +} diff --git a/src/olbProcessors/saturatedFluxPostProcessor2D.hh b/src/olbProcessors/saturatedFluxPostProcessor2D.hh new file mode 100644 index 0000000..47e46a5 --- /dev/null +++ b/src/olbProcessors/saturatedFluxPostProcessor2D.hh @@ -0,0 +1,74 @@ +#include "saturatedFluxPostProcessor2D.h" + +#include +#include + +namespace olb { + +//////// SaturationFluxProcessor2D //////////////////////////// + +template +SaturatedFluxProcessor2D:: +SaturatedFluxProcessor2D(int x0_, int x1_, int y0_, int y1_, int discreteNormalX_, int discreteNormalY_, T *Vmax_, T *Km_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), + discreteNormalX(discreteNormalX_), discreteNormalY(discreteNormalY_), + Vmax(Vmax_), Km(Km_) +{ + this->getName() = "SaturatedFluxProcessor2D"; +} + +template +void SaturatedFluxProcessor2D:: +processSubDomain(BlockLattice& blockLattice, int x0_, int x1_, int y0_, int y1_) +{ + int newX0, newX1, newY0, newY1; + if ( util::intersect ( + x0, x1, y0, y1, + x0_, x1_, y0_, y1_, + newX0, newX1, newY0, newY1 ) ) { + + T u[2]; + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + T tempFluid = blockLattice.get(iX-discreteNormalX, iY-discreteNormalY).computeRho(); + T flux = ((*Vmax) * tempFluid) / ((*Km) + tempFluid); + u[0] = flux*discreteNormalX; + u[1] = flux*discreteNormalY; + blockLattice.get(iX,iY).defineU(u); + } + } + } +} + +template +void SaturatedFluxProcessor2D:: +process(BlockLattice& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1); +} + +template +SaturatedFluxProcessorGenerator2D:: +SaturatedFluxProcessorGenerator2D(int x0_, int x1_, int y0_, int y1_, int discreteNormalX_, + int discreteNormalY_, T *Vmax_, T *Km_) + : PostProcessorGenerator2D(x0_, x1_, y0_, y1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), Vmax(Vmax_), Km(Km_) +{ } + +template +PostProcessor2D* +SaturatedFluxProcessorGenerator2D::generate() const +{ + return new SaturatedFluxProcessor2D + (this->x0, this->x1, this->y0, this->y1, discreteNormalX, discreteNormalY, Vmax, Km); +} + +template +PostProcessorGenerator2D* +SaturatedFluxProcessorGenerator2D::clone() const +{ + return new SaturatedFluxProcessorGenerator2D + (this->x0, this->x1, this->y0, this->y1, discreteNormalX, discreteNormalY, Vmax, Km); +} + +} diff --git a/src/olbProcessors/setFunctionalRegularizedHeatFlux.h b/src/olbProcessors/setFunctionalRegularizedHeatFlux.h new file mode 100644 index 0000000..fd4da26 --- /dev/null +++ b/src/olbProcessors/setFunctionalRegularizedHeatFlux.h @@ -0,0 +1,22 @@ +#pragma once + +#include +#include + +namespace olb { + +///Initialising the FunctionalRegularizedHeatFluxBoundary on the superLattice domain +///This is an advection diffusion boundary -->MixinDynamics = AdvectionDiffusionRLBdynamics +template> +void setFunctionalRegularizedHeatFluxBoundary(SuperLattice& sLattice,T omega, SuperGeometry& superGeometry, int material, T *Vmax=nullptr, T *Km=nullptr); + +///Initialising the RegularizedHeatFluxBoundary on the superLattice domain +template +void setFunctionalRegularizedHeatFluxBoundary(SuperLattice& sLattice,T omega, FunctorPtr>&& indicator, T *Vmax=nullptr, T *Km=nullptr); + + +///Set RegularizedHeatFluxBoundary for indicated cells inside the block domain +template +void setFunctionalRegularizedHeatFluxBoundary(BlockLattice& block, T omega, BlockIndicatorF2D& indicator, T *Vmax, T *Km); + +}//namespace olb diff --git a/src/olbProcessors/setFunctionalRegularizedHeatFlux.hh b/src/olbProcessors/setFunctionalRegularizedHeatFlux.hh new file mode 100644 index 0000000..7c93bd9 --- /dev/null +++ b/src/olbProcessors/setFunctionalRegularizedHeatFlux.hh @@ -0,0 +1,86 @@ +#include "setFunctionalRegularizedHeatFlux.h" + +#include +#include + +namespace olb { + +/** +* Class for setting a regularized head flux boundary condition on +* advection diffusion fields, depending on the temperature (concentration) +* value of the neighbouring node(s). +*/ + +// Initialising the FunctionalRegularizedHeatFluxBoundary on the superLattice domain +template +void setFunctionalRegularizedHeatFluxBoundary(SuperLattice& sLattice, T omega, SuperGeometry& superGeometry, int material, T *Vmax, T *Km) +{ + setFunctionalRegularizedHeatFluxBoundary(sLattice, omega, superGeometry.getMaterialIndicator(material), Vmax, Km); +} + +// Initializing the FunctionalRegularizedHeatFluxBoundary on the superLattice domain +template +void setFunctionalRegularizedHeatFluxBoundary(SuperLattice& sLattice, T omega, FunctorPtr>&& indicator, T *Vmax, T *Km) +{ + int _overlap = 1; + for (int iCloc = 0; iCloc < sLattice.getLoadBalancer().size(); ++iCloc) { + setFunctionalRegularizedHeatFluxBoundary(sLattice.getBlock(iCloc), omega, + indicator->getBlockIndicatorF(iCloc), Vmax, Km); + } + // Adds needed Cells to the Communicator _commBC in SuperLattice + addPoints2CommBC(sLattice, std::forward(indicator), _overlap); +} + + +// Set FunctionalRegularizedHeatFluxBoundary for indicated cells inside the block domain +template +void setFunctionalRegularizedHeatFluxBoundary(BlockLattice& block, T omega, BlockIndicatorF2D& indicator, T *Vmax, T *Km) +{ + using namespace boundaryhelper; + auto& blockGeometryStructure = indicator.getBlockGeometry(); + OstreamManager clout(std::cout, "setFunctionalRegularizedHeatFluxBoundary"); + /* + *x0,x1,y0,y1 Range of cells to be traversed + **/ + int x0 = 0; + int y0 = 0; + int x1 = blockGeometryStructure.getNx()-1; + int y1 = blockGeometryStructure.getNy()-1; + std::vector discreteNormal(3,0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + Dynamics* dynamics = nullptr; + if (indicator(iX, iY)) { + discreteNormal = indicator.getBlockGeometry().getStatistics().getType(iX,iY); + if (discreteNormal[0] == 0) { // set momenta, dynamics and posprocessor on indicated FunctionalRegularizedHeatFluxBoundary cells + dynamics = block.getDynamics(PlainMixinDynamicsForDirectionOrientationMomenta::construct(Vector(discreteNormal.data() + 1))); + setBoundary(block, iX, iY, dynamics); + + auto fluxPostProcessor = std::unique_ptr>{ + new SaturatedFluxProcessorGenerator2D( + iX, iX, iY, iY, discreteNormal[1], discreteNormal[2], Vmax, Km ) + }; + + block.addPostProcessor(*fluxPostProcessor); + } + else if (discreteNormal[0] == 1) { // set momenta, dynamics and postProcessor on indicated FunctionalRegularizedHeatFluxBoundary external corner cells + dynamics = block.getDynamics(NormalMixinDynamicsForPlainMomenta::construct(Vector(discreteNormal.data() + 1))); + setBoundary(block, iX, iY, dynamics); + } + if (dynamics) { + dynamics->getParameters(block).template set(omega); + } + } + } + } +} + + + +} diff --git a/src/porting/jsonPorter.hh b/src/porting/jsonPorter.hh index c7f4bfc..1f7dc25 100644 --- a/src/porting/jsonPorter.hh +++ b/src/porting/jsonPorter.hh @@ -113,10 +113,23 @@ void simulationFromJSON(json jsonString, arch::Network* network_, sim::Simula if (platform == sim::Platform::Continuous) { readSimulators(jsonString, simulation, network_); network_->sortGroups(); + } else if (platform == sim::Platform::Mixing) { + readMixingModel(jsonString, simulation); + readSpecies(jsonString, simulation); + readMixtures(jsonString, simulation); + readMixtureInjections(jsonString, simulation, activeFixture); + readSimulators(jsonString, simulation, network_); + network_->sortGroups(); + } else if (platform == sim::Platform::Ooc) { + readMixingModel(jsonString, simulation); + readSpecies(jsonString, simulation); + readMixtures(jsonString, simulation); + readMixtureInjections(jsonString, simulation, activeFixture); + readTissues(jsonString, simulation); + readSimulators(jsonString, simulation, network_); + network_->sortGroups(); } else if (platform == sim::Platform::BigDroplet) { throw std::invalid_argument("Droplet simulations are currently only supported for Abstract simulations."); - } else if (platform == sim::Platform::Mixing) { - throw std::invalid_argument("Mixing simulations are currently only supported for Abstract simulations."); } else { throw std::invalid_argument("Invalid platform for Hybrid simulation. Please select one of the following:\n\tcontinuous"); } diff --git a/src/porting/jsonReaders.h b/src/porting/jsonReaders.h index 914cdbe..4d00b14 100644 --- a/src/porting/jsonReaders.h +++ b/src/porting/jsonReaders.h @@ -102,6 +102,9 @@ void readDroplets (json jsonString, sim::Simulation& simulation); template void readSpecies (json jsonString, sim::Simulation& simulation); +template +void readTissues (json jsonString, sim::Simulation& simulation); + template void readMixtures (json jsonString, sim::Simulation& simulation); diff --git a/src/porting/jsonReaders.hh b/src/porting/jsonReaders.hh index b20cccc..d24e32d 100644 --- a/src/porting/jsonReaders.hh +++ b/src/porting/jsonReaders.hh @@ -70,7 +70,7 @@ template sim::Platform readPlatform(json jsonString, sim::Simulation& simulation) { sim::Platform platform = sim::Platform::Continuous; if (!jsonString["simulation"].contains("platform")) { - throw std::invalid_argument("Please define a platform. The following platforms are possible:\nContinuous\nBigDroplet\nMixing"); + throw std::invalid_argument("Please define a platform. The following platforms are possible:\nContinuous\nBigDroplet\nMixing\nOoc"); } if (jsonString["simulation"]["platform"] == "Continuous") { platform = sim::Platform::Continuous; @@ -78,8 +78,10 @@ sim::Platform readPlatform(json jsonString, sim::Simulation& simulation) { platform = sim::Platform::BigDroplet; } else if (jsonString["simulation"]["platform"] == "Mixing") { platform = sim::Platform::Mixing; + } else if (jsonString["simulation"]["platform"] == "Ooc") { + platform = sim::Platform::Ooc; } else { - throw std::invalid_argument("Platform is invalid. The following platforms are possible:\nContinuous\nBigDroplet\nMixing"); + throw std::invalid_argument("Platform is invalid. The following platforms are possible:\nContinuous\nBigDroplet\nMixing\nOoc"); } simulation.setPlatform(platform); return platform; @@ -149,6 +151,32 @@ void readSpecies(json jsonString, sim::Simulation& simulation) { } } +template +void readTissues(json jsonString, sim::Simulation& simulation) { + for (auto& tissue : jsonString["simulation"]["tissues"]) { + if (tissue.contains("species") && tissue.contains("Vmax") && tissue.contains("kM")) { + if (tissue["species"].size() == tissue["Vmax"].size() && tissue["species"].size() == tissue["kM"].size()) { + int counter = 0; + std::unordered_map*> species; + std::unordered_map Vmax; + std::unordered_map kM; + for (auto& specieId : tissue["species"]) { + auto specie_ptr = simulation.getSpecie(specieId); + species.try_emplace(specieId, specie_ptr); + Vmax.try_emplace(specieId, tissue["Vmax"][counter]); + kM.try_emplace(specieId, tissue["kM"][counter]); + counter++; + } + simulation.addTissue(species, Vmax, kM); + } else { + throw std::invalid_argument("Please define a Vmax and kM value for each species."); + } + } else { + throw std::invalid_argument("Wrongly defined specie. Please provide following information for species:\ndiffusivity\nsaturationConcentration"); + } + } +} + template void readMixtures(json jsonString, sim::Simulation& simulation) { for (auto& mixture : jsonString["simulation"]["mixtures"]) { @@ -245,6 +273,28 @@ void readSimulators(json jsonString, sim::Simulation& simulation, arch::Netwo charPhysVelocity, alpha, resolution, epsilon, tau); simulator->setVtkFolder(vtkFolder); } + else if (simulator["Type"] == "Mixing") + { + std::unordered_map*> species; + for (auto& [specieId, speciePtr] : simulation.getSpecies()) { + species.try_emplace(specieId, speciePtr.get()); + } + auto simulator = simulation.addLbmMixingSimulator(name, stlFile, network->getModule(moduleId), species, + Openings, charPhysLength, charPhysVelocity, alpha, resolution, epsilon, tau); + simulator->setVtkFolder(vtkFolder); + } + else if (simulator["Type"] == "Organ") + { + std::unordered_map*> species; + for (auto& [specieId, speciePtr] : simulation.getSpecies()) { + species.try_emplace(specieId, speciePtr.get()); + } + std::string organStlFile = simulator["organStlFile"]; + int tissueId = simulator["tissue"]; + auto simulator = simulation.addLbmOocSimulator(name, stlFile, tissueId, organStlFile, network->getModule(moduleId), species, + Openings, charPhysLength, charPhysVelocity, alpha, resolution, epsilon, tau); + simulator->setVtkFolder(vtkFolder); + } else if(simulator["Type"] == "ESS_LBM") { #ifdef USE_ESSLBM diff --git a/src/simulation/CFDSim.h b/src/simulation/CFDSim.h index 02b6a71..54922f1 100644 --- a/src/simulation/CFDSim.h +++ b/src/simulation/CFDSim.h @@ -20,4 +20,10 @@ class CFDSimulator; template bool conductCFDSimulation(const std::unordered_map>>& cfdSimulators, int iteration); + template + void coupleNsAdLattices(const std::unordered_map>>& cfdSimulators); + + template + bool conductADSimulation(const std::unordered_map>>& cfdSimulators); + } // namespace sim \ No newline at end of file diff --git a/src/simulation/CFDSim.hh b/src/simulation/CFDSim.hh index bca12dc..d855f32 100644 --- a/src/simulation/CFDSim.hh +++ b/src/simulation/CFDSim.hh @@ -27,4 +27,49 @@ namespace sim { return allConverge; } + template + void coupleNsAdLattices(const std::unordered_map>>& cfdSimulators) { + + // loop through modules and perform the coupling between the NS and AD lattices + for (const auto& cfdSimulator : cfdSimulators) { + + // Assertion that the current module is of lbm type, and can conduct CFD simulations. + #ifndef USE_ESSLBM + assert(cfdSimulator.second->getModule()->getModuleType() == arch::ModuleType::LBM); + #elif USE_ESSLBM + assert(cfdSimulator.second->getModule()->getModuleType() == arch::ModuleType::ESS_LBM); + throw std::runtime_error("Coupling between NS and AD fields not defined for ESS LBM."); + #endif + + cfdSimulator.second->executeCoupling(); + + } + } + + template + bool conductADSimulation(const std::unordered_map>>& cfdSimulators) { + + bool allConverge = true; + + // loop through modules and perform the collide and stream operations + for (const auto& cfdSimulator : cfdSimulators) { + + // Assertion that the current module is of lbm type, and can conduct CFD simulations. + #ifndef USE_ESSLBM + assert(cfdSimulator.second->getModule()->getModuleType() == arch::ModuleType::LBM); + #elif USE_ESSLBM + assert(cfdSimulator.second->getModule()->getModuleType() == arch::ModuleType::ESS_LBM); + throw std::runtime_error("Simulation of Advection Diffusion not defined for ESS LBM."); + #endif + cfdSimulator.second->adSolve(); + + if (!cfdSimulator.second->hasAdConverged()) { + allConverge = false; + } + + } + + return allConverge; + } + } // namespace sim diff --git a/src/simulation/CMakeLists.txt b/src/simulation/CMakeLists.txt index d8fa674..f23e07e 100644 --- a/src/simulation/CMakeLists.txt +++ b/src/simulation/CMakeLists.txt @@ -9,6 +9,7 @@ set(SOURCE_LIST ResistanceModels.hh Simulation.hh Specie.hh + Tissue.hh ) set(HEADER_LIST @@ -22,6 +23,7 @@ set(HEADER_LIST ResistanceModels.h Simulation.h Specie.h + Tissue.h ) target_sources(${TARGET_NAME} PUBLIC ${SOURCE_LIST} ${HEADER_LIST}) diff --git a/src/simulation/MembraneModels.h b/src/simulation/MembraneModels.h new file mode 100644 index 0000000..c22c032 --- /dev/null +++ b/src/simulation/MembraneModels.h @@ -0,0 +1,61 @@ +/** + * @file MembraneModels.h + */ + +#pragma once + +#include + +namespace sim { + +/** + * @brief Virtual class that describes the basic functionality for mixing models. +*/ +template +class MembraneModel { +protected: + +public: + + /** + * @brief Constructor of a membrane model. + */ + MembraneModel(); + +}; + +/** + * @brief Class that describes the poreResistance membrane model. This membrane model derives the permeability from pore geometry. +*/ +template +class PermeabilityMembraneModel : public MembraneModel { + +private: + +public: + + /** + * @brief Constructor of a poreResistance membrane model. + */ + PermeabilityMembraneModel(); + +}; + +/** + * @brief Class that describes the poreResistance membrane model. This membrane model derives the permeability from pore geometry. +*/ +template +class PoreResistanceMembraneModel : public MembraneModel { + +private: + +public: + + /** + * @brief Constructor of a poreResistance membrane model. + */ + PoreResistanceMembraneModel(); + +}; + +} // namespace sim \ No newline at end of file diff --git a/src/simulation/MembraneModels.hh b/src/simulation/MembraneModels.hh new file mode 100644 index 0000000..5ddac64 --- /dev/null +++ b/src/simulation/MembraneModels.hh @@ -0,0 +1,26 @@ +#include "MembraneModels.h" + +#include + +namespace sim { + +template +MembraneModel::MembraneModel() { } + +//========================================================================================= +//============================= permeability Membrane =================================== +//========================================================================================= + +template +PermeabilityMembraneModel::PermeabilityMembraneModel() : MembraneModel() { } + + +//========================================================================================= +//============================ poreResistance Membrane ================================== +//========================================================================================= + +template +PoreResistanceMembraneModel::PoreResistanceMembraneModel() : MembraneModel() { } + + +} /// namespace sim diff --git a/src/simulation/MixingModels.h b/src/simulation/MixingModels.h index 1dd39b0..7fbd0c3 100644 --- a/src/simulation/MixingModels.h +++ b/src/simulation/MixingModels.h @@ -89,6 +89,11 @@ class MixingModel { */ MixingModel(); + /** + * @brief Propagate all the species through a network for a steady-state simulation + */ + virtual void propagateSpecies(arch::Network* network, Simulation* sim) = 0; + /** * @brief Returns the current minimal timestep. * @return The minimal timestep in s. @@ -157,7 +162,7 @@ class InstantaneousMixingModel : public MixingModel { private: std::unordered_map>> mixtureInflowAtNode; ///< Unordered map to track mixtures flowing into nodes > - std::unordered_map mixtureOutflowAtNode; ///< Unordered map to track mixtures flowing out of nodes. + std::unordered_map mixtureOutflowAtNode; ///< Unordered map to track mixtures flowing out of nodes . std::unordered_map totalInflowVolumeAtNode; ///< Unordered map to track the total volumetric flow entering a node. std::unordered_map createMixture; ///< Unordered map to track whether a new mixture is created at a node. @@ -207,6 +212,30 @@ class InstantaneousMixingModel : public MixingModel { */ void clean(arch::Network* network); + /** + * @brief Propagate all the species through a network for a steady-state simulation + */ + void propagateSpecies(arch::Network* network, Simulation* sim) override; + + /** + * @brief From the mixtureInjections and CFD simulators, generate temporary mxtures that + * flow into the network at correspondingnode entry points. + * @param[in] sim Pointer to the simulation. + */ + void initNodeOutflow(Simulation* sim, std::vector>& tmpMixtures); + + /** + * @brief Propagate the mixtures through the corresponding channel entirely, without considering time steps + */ + void channelPropagation(arch::Network* network); + + /** + * @brief From the node's inflows, generate the node outflow + */ + bool updateNodeOutflow(Simulation* sim, std::vector>& tmpMixtures); + + void storeConcentrations(Simulation* sim, const std::vector>& tmpMixtures); + /** * @brief Print all mixtures and their positions in the network. */ @@ -253,6 +282,11 @@ class DiffusionMixingModel : public MixingModel { void generateInflows(T timeStep, arch::Network* network, Simulation* sim, std::unordered_map>>& mixtures); void topologyAnalysis(arch::Network* network, int nodeId); + + /** + * @brief Propagate all the species through a network for a steady-state simulation + */ + void propagateSpecies(arch::Network* network, Simulation* sim) override; void printTopology(); diff --git a/src/simulation/MixingModels.hh b/src/simulation/MixingModels.hh index 0ee3388..f2fd294 100644 --- a/src/simulation/MixingModels.hh +++ b/src/simulation/MixingModels.hh @@ -61,7 +61,168 @@ void MixingModel::injectMixtureInEdge(int mixtureId, int channelId) { } template -InstantaneousMixingModel::InstantaneousMixingModel() { } +InstantaneousMixingModel::InstantaneousMixingModel() : MixingModel() { } + +template +void InstantaneousMixingModel::propagateSpecies(arch::Network* network, Simulation* sim) { + + std::vector> tmpMixtures; + + // Define total inflow volume at nodes + for (auto& [nodeId, node] : network->getNodes()) { + for (auto& channel : network->getChannelsAtNode(nodeId) ) { + // Check if the channel flows into the node + if ((channel->getFlowRate() > 0.0 && channel->getNodeB() == nodeId) || (channel->getFlowRate() < 0.0 && channel->getNodeA() == nodeId)) { + T inflowVolume = std::abs(channel->getFlowRate()); + auto [iterator, inserted] = totalInflowVolumeAtNode.try_emplace(nodeId, inflowVolume); + if (!inserted) { + iterator->second += inflowVolume; + } + } + } + } + + // Initial node outflow from mixtureInjections and CFD simulators, stored in mixtureOutflowAtNode + initNodeOutflow(sim, tmpMixtures); + + // Propagate the mixtures through the entire channel, without considering time steps + channelPropagation(network); + + bool inflowUpdated = true; + while (inflowUpdated) { + // From node inflow, generate the node's outflow + inflowUpdated = updateNodeOutflow(sim, tmpMixtures); + // Propagate the mixtures through the entire channel + channelPropagation(network); + } + + // Store the concentrations of the final state in the concentration buffer of olbMixingSolver. + storeConcentrations(sim, tmpMixtures); + + clean(network); + +} + +template +void InstantaneousMixingModel::initNodeOutflow(Simulation* sim, std::vector>& tmpMixtures) { + // Add mixture injections + for (auto& [key, mixtureInjection] : sim->getMixtureInjections()) { + int tmpMixtureIndex = tmpMixtures.size(); + int nodeId = mixtureInjection->getInjectionChannel()->getFlowRate() >= 0.0 ? + mixtureInjection->getInjectionChannel()->getNodeB() : mixtureInjection->getInjectionChannel()->getNodeA(); + tmpMixtures.push_back(Mixture(*sim->getMixture(mixtureInjection->getMixtureId()))); + mixtureOutflowAtNode.try_emplace(nodeId, tmpMixtureIndex); + } + // Add CFD Simulator outflows + for (auto& [key, cfdSimulator] : sim->getCFDSimulators()) { + for (auto& [nodeId, opening] : cfdSimulator->getOpenings()) { + // If the node is an outflow + if (cfdSimulator->getFlowRates().at(nodeId) < 0.0) { + int tmpMixtureIndex = tmpMixtures.size(); + int tmpMixtureId = std::numeric_limits::max(); + std::unordered_map*> species; + std::unordered_map speciesConcentrations(cfdSimulator->getConcentrations().at(nodeId)); + for (auto& [speciesId, concentration] : cfdSimulator->getConcentrations().at(nodeId)) { + species.try_emplace(speciesId, sim->getSpecie(speciesId)); + } + Mixture tmpMixture = Mixture(tmpMixtureId, species, speciesConcentrations, sim->getContinuousPhase()); + tmpMixtures.push_back(tmpMixture); + mixtureOutflowAtNode.try_emplace(nodeId, tmpMixtureIndex); + } + } + } +} + +template +void InstantaneousMixingModel::channelPropagation(arch::Network* network) { + for (auto& [nodeId, mixtureId] : mixtureOutflowAtNode) { + for (auto& channel : network->getChannelsAtNode(nodeId)) { + // Find the nodeId that is across the channel + int oppositeNode; + if (channel->getFlowRate() > 0.0 && channel->getNodeA() == nodeId) { + oppositeNode = channel->getNodeB(); + } else if (channel->getFlowRate() < 0.0 && channel->getNodeB() == nodeId) { + oppositeNode = channel->getNodeA(); + } else { + continue; + } + // Update the mixture inflow at the node across the channel + MixtureInFlow mixtureInflow = {mixtureId, std::abs(channel->getFlowRate())}; + auto [iterator, inserted] = mixtureInflowAtNode.try_emplace(oppositeNode, std::vector>(1, mixtureInflow)); + if (!inserted) { + mixtureInflowAtNode.at(oppositeNode).push_back(mixtureInflow); + } + } + } +} + +template +bool InstantaneousMixingModel::updateNodeOutflow(Simulation* sim, std::vector>& tmpMixtures) { + bool updated = false; + for (auto& [nodeId, mixtureInflowList] : mixtureInflowAtNode) { + bool createMixture = false; + std::unordered_map*> speciePtrs; + std::unordered_map newConcentrations; + for (auto& mixtureInflow : mixtureInflowList) { + for (auto& [specieId, oldConcentration] : tmpMixtures[mixtureInflow.mixtureId].getSpecieConcentrations()) { + speciePtrs.try_emplace(specieId, sim->getSpecie(specieId)); + T newConcentration = oldConcentration * mixtureInflow.inflowVolume / totalInflowVolumeAtNode.at(nodeId); + auto [iterator, inserted] = newConcentrations.try_emplace(specieId, newConcentration); + if (!inserted) { + iterator->second = iterator->second + newConcentration; + } + } + if (mixtureInflow.mixtureId != mixtureInflowList[0].mixtureId) { + createMixture = true; + } + } + int outflowMixtureId; + Mixture newMixture (tmpMixtures.size(), speciePtrs, newConcentrations, sim->getContinuousPhase()); + if (createMixture) { + outflowMixtureId = tmpMixtures.size(); + } else { + outflowMixtureId = mixtureInflowList[0].mixtureId; + } + if (!mixtureOutflowAtNode.count(nodeId)) { + if (createMixture) { + tmpMixtures.push_back(newMixture); + } + mixtureOutflowAtNode.try_emplace(nodeId, outflowMixtureId); + updated = true; + + } else { + // Check if the ouflow gets updated or not + if (tmpMixtures[mixtureOutflowAtNode.at(nodeId)] == tmpMixtures[outflowMixtureId]) { + continue; + } else { + if (createMixture) { + tmpMixtures.push_back(newMixture); + } + mixtureOutflowAtNode.at(nodeId) = outflowMixtureId; + updated = true; + } + } + } + return updated; +} + +template +void InstantaneousMixingModel::storeConcentrations(Simulation* sim, const std::vector>& tmpMixtures) { + for (auto& [key, cfdSimulator] : sim->getCFDSimulators()) { + std::unordered_map> concentrations = cfdSimulator->getConcentrations(); + for (auto& [nodeId, opening] : cfdSimulator->getOpenings()) { + // If the node is an inflow + if (cfdSimulator->getFlowRates().at(nodeId) > 0.0) { + if (mixtureOutflowAtNode.count(nodeId)) { + for (auto& [specieId, specieConcentration] : tmpMixtures[mixtureOutflowAtNode.at(nodeId)].getSpecieConcentrations()) { + concentrations.at(nodeId).at(specieId) = specieConcentration; + } + } + } + } + cfdSimulator->storeConcentrations(concentrations); + } +} template void InstantaneousMixingModel::updateMixtures(T timeStep, arch::Network* network, Simulation* sim, std::unordered_map>>& mixtures) { @@ -502,7 +663,12 @@ void DiffusionMixingModel::topologyAnalysis( arch::Network* network, int n } template -std::tuple, std::vector, T>DiffusionMixingModel::getAnalyticalSolutionConstant(T channelLength, T channelWidth, int resolution, T pecletNr, const std::vector>& parameters) { +void DiffusionMixingModel::propagateSpecies(arch::Network* network, Simulation* sim) { + // TODO +} + +template +std::tuple, std::vector, T> DiffusionMixingModel::getAnalyticalSolutionConstant(T channelLength, T channelWidth, int resolution, T pecletNr, const std::vector>& parameters) { T a_0 = 0.0; std::vector segmentedResult; diff --git a/src/simulation/Mixture.h b/src/simulation/Mixture.h index 34920eb..a9613c0 100644 --- a/src/simulation/Mixture.h +++ b/src/simulation/Mixture.h @@ -90,6 +90,11 @@ class Mixture { */ Mixture(int id, std::unordered_map*> species, std::unordered_map specieConcentrations, Fluid* carrierFluid); + /** + * @brief Overload == operator to check if mixtures are equal or not + */ + bool operator== (const Mixture &t); + /** * @brief Get the id of this mixture * @return Unique identifier of the mixture. diff --git a/src/simulation/Mixture.hh b/src/simulation/Mixture.hh index cc8a145..3cb66c7 100644 --- a/src/simulation/Mixture.hh +++ b/src/simulation/Mixture.hh @@ -21,6 +21,17 @@ Mixture::Mixture(int id, std::unordered_map*> species, std::un viscosity(carrierFluid->getViscosity()), density(carrierFluid->getDensity()), largestMolecularSize(0.0) { } +template +bool Mixture::operator== (const Mixture &t) { + if (species == t.species && + specieConcentrations == t.specieConcentrations && + viscosity == t.viscosity && + density == t.density) { + return true; + } + return false; +} + template int Mixture::getId() const { return id; diff --git a/src/simulation/Simulation.h b/src/simulation/Simulation.h index dca9fee..c4a2b62 100644 --- a/src/simulation/Simulation.h +++ b/src/simulation/Simulation.h @@ -67,9 +67,18 @@ class Fluid; template class lbmSimulator; +template +class lbmMixingSimulator; + +template +class lbmOocSimulator; + template class essLbmSimulator; +template +class MembraneModel; + template class MixingModel; @@ -85,6 +94,9 @@ class ResistanceModel; template class Specie; +template +class Tissue; + enum class Type { Abstract, ///< A simulation in the 1D abstraction level Hybrid, ///< A simulation combining the 1D and CFD abstraction levels @@ -94,7 +106,8 @@ enum class Type { enum class Platform { Continuous, ///< A simulation with a single continuous fluid. BigDroplet, ///< A simulation with droplets filling a channel cross-section - Mixing ///< A simulation wit multiple miscible fluids. + Mixing, ///< A simulation with multiple miscible fluids. + Ooc ///< A simulation with organic tissue }; /** @@ -111,11 +124,13 @@ class Simulation { std::unordered_map>> fluids; ///< Fluids specified for the simulation. std::unordered_map>> droplets; ///< Droplets which are simulated in droplet simulation. std::unordered_map>> species; ///< Species specified for the simulation. + std::unordered_map>> tissues; ///< Tissues specified for the simulation. std::unordered_map>> dropletInjections; ///< Injections of droplets that should take place during a droplet simulation. std::unordered_map>> mixtures; ///< Mixtures present in the simulation. std::unordered_map>> mixtureInjections; ///< Injections of fluids that should take place during the simulation. std::unordered_map>> cfdSimulators; ResistanceModel* resistanceModel; ///< The resistance model used for the simulation. + MembraneModel* membraneModel; ///< The membrane model used for an OoC simulation. MixingModel* mixingModel; ///< The mixing model used for a mixing simulation. int continuousPhase = 0; ///< Fluid of the continuous phase. int iteration = 0; @@ -208,6 +223,15 @@ class Simulation { */ Specie* addSpecie(T diffusivity, T satConc); + /** + * @brief Create tissue. + * @param[in] species Map of Species that interacts with this tissue. + * @param[in] Vmax + * @param[in] kM + * @return Pointer to created tissue. + */ + Tissue* addTissue(std::unordered_map*> species, std::unordered_map Vmax, std::unordered_map kM); + /** * @brief Create mixture. * @param[in] specieConcentrations unordered map of specie id and corresponding concentration. @@ -289,6 +313,45 @@ class Simulation { lbmSimulator* addLbmSimulator(std::string name, std::string stlFile, std::shared_ptr> module, std::unordered_map> openings, T charPhysLength, T charPhysVelocity, T alpha, T resolution, T epsilon, T tau); + /** + * @brief Adds a new module to the network. + * @param[in] name Name of the module. + * @param[in] stlFile Location of the stl file that gives the geometry of the domain. + * @param[in] module Shared pointer to the module on which this solver acts. + * @param[in] species Map of specieIds and speciePtrs of the species simulated in the AD fields of this simulator. + * @param[in] openings Map of openings corresponding to the nodes. + * @param[in] charPhysLength Characteristic physical length of this simulator. + * @param[in] charPhysVelocity Characteristic physical velocity of this simulator. + * @param[in] alpha Relaxation parameter for this simulator. + * @param[in] resolution Resolution of this simulator. + * @param[in] epsilon Error tolerance for convergence criterion of this simulator. + * @param[in] tau Relaxation time of this simulator (0.5 < tau < 2.0). + * @return Pointer to the newly created module. + */ + lbmMixingSimulator* addLbmMixingSimulator(std::string name, std::string stlFile, std::shared_ptr> module, std::unordered_map*> species, + std::unordered_map> openings, T charPhysLength, T charPhysVelocity, T alpha, T resolution, T epsilon, T tau); + + + /** + * @brief Adds a new module to the network. + * @param[in] name Name of the module. + * @param[in] stlFile Location of the stl file that gives the geometry of the domain. + * @param[in] tissueId The Id of the tissue that the organ in this nodule consists of. + * @param[in] organStlFile The location of the stl file describing the geometry of the organ. + * @param[in] module Shared pointer to the module on which this solver acts. + * @param[in] species Map of specieIds and speciePtrs of the species simulated in the AD fields of this simulator. + * @param[in] openings Map of openings corresponding to the nodes. + * @param[in] charPhysLength Characteristic physical length of this simulator. + * @param[in] charPhysVelocity Characteristic physical velocity of this simulator. + * @param[in] alpha Relaxation parameter for this simulator. + * @param[in] resolution Resolution of this simulator. + * @param[in] epsilon Error tolerance for convergence criterion of this simulator. + * @param[in] tau Relaxation time of this simulator (0.5 < tau < 2.0). + * @return Pointer to the newly created module. + */ + lbmOocSimulator* addLbmOocSimulator(std::string name, std::string stlFile, int tissueId, std::string organStlFile, std::shared_ptr> module, std::unordered_map*> species, + std::unordered_map> openings, T charPhysLength, T charPhysVelocity, T alpha, T resolution, T epsilon, T tau); + /** * @brief Adds a new module to the network. * @param[in] name Name of the module. @@ -431,6 +494,25 @@ class Simulation { */ MixtureInjection* getMixtureInjection(int injectionId); + /** + * @brief Get injection + * @return Reference to the unordered map of MixtureInjections + */ + std::unordered_map>>& getMixtureInjections(); + + /** + * @brief Get injection + * @param simulatorId The id of the injection + * @return Pointer to injection with the corresponding id. + */ + CFDSimulator* getCFDSimulator(int simulatorId); + + /** + * @brief Get injection + * @return Reference to the unordered map of MixtureInjections + */ + std::unordered_map>>& getCFDSimulators(); + /** * @brief Get the continuous phase. * @return Fluid if the continuous phase or nullptr if no continuous phase is specified. diff --git a/src/simulation/Simulation.hh b/src/simulation/Simulation.hh index 65a21cc..8fb8553 100644 --- a/src/simulation/Simulation.hh +++ b/src/simulation/Simulation.hh @@ -35,6 +35,15 @@ namespace sim { return result.first->second.get(); } + template + Tissue* Simulation::addTissue(std::unordered_map*> species, std::unordered_map Vmax, std::unordered_map kM) { + auto id = tissues.size(); + + auto result = tissues.insert_or_assign(id, std::make_shared>(id, species, Vmax, kM)); + + return result.first->second.get(); + } + template Mixture* Simulation::addMixture(std::unordered_map specieConcentrations) { auto id = mixtures.size(); @@ -207,6 +216,43 @@ namespace sim { } } + template + lbmMixingSimulator* Simulation::addLbmMixingSimulator(std::string name, std::string stlFile, std::shared_ptr> module, std::unordered_map*> species, + std::unordered_map> openings, T charPhysLength, T charPhysVelocity, T alpha, T resolution, T epsilon, T tau) + { + std::cout << "Trying to add a mixing simulator" << std::endl; + if (resistanceModel != nullptr) { + // create Simulator + auto id = cfdSimulators.size(); + auto addCfdSimulator = new lbmMixingSimulator(id, name, stlFile, module, species, openings, resistanceModel, charPhysLength, charPhysVelocity, alpha, resolution, epsilon, tau); + + // add Simulator + cfdSimulators.try_emplace(id, addCfdSimulator); + + return addCfdSimulator; + } else { + throw std::invalid_argument("Attempt to add CFD Simulator without valid resistanceModel."); + } + } + + template + lbmOocSimulator* Simulation::addLbmOocSimulator(std::string name, std::string stlFile, int tissueId, std::string organStlFile, std::shared_ptr> module, std::unordered_map*> species, + std::unordered_map> openings, T charPhysLength, T charPhysVelocity, T alpha, T resolution, T epsilon, T tau) + { + if (resistanceModel != nullptr) { + // create Simulator + auto id = cfdSimulators.size(); + auto addCfdSimulator = new lbmOocSimulator(id, name, stlFile, tissues.at(tissueId), organStlFile, module, species, openings, resistanceModel, charPhysLength, charPhysVelocity, alpha, resolution, epsilon, tau); + + // add Simulator + cfdSimulators.try_emplace(id, addCfdSimulator); + + return addCfdSimulator; + } else { + throw std::invalid_argument("Attempt to add CFD Simulator without valid resistanceModel."); + } + } + template essLbmSimulator* Simulation::addEssLbmSimulator(std::string name, std::string stlFile, std::shared_ptr> module, std::unordered_map> openings, T charPhysLength, T charPhysVelocity, T alpha, T resolution, T epsilon, T tau) @@ -361,6 +407,21 @@ namespace sim { return mixtureInjections.at(injectionId).get(); } + template + std::unordered_map>>& Simulation::getMixtureInjections() { + return mixtureInjections; + } + + template + CFDSimulator* Simulation::getCFDSimulator(int simulatorId) { + return cfdSimulators.at(simulatorId).get(); + } + + template + std::unordered_map>>& Simulation::getCFDSimulators() { + return cfdSimulators; + } + template Fluid* Simulation::getContinuousPhase() { return fluids[continuousPhase].get(); @@ -484,35 +545,109 @@ namespace sim { // Continuous Hybrid simulation if (this->simType == Type::Hybrid && this->platform == Platform::Continuous) { - if (network->getModules().size() > 0 ) { - bool allConverged = false; - bool pressureConverged = false; + + // Catch runtime error, not enough CFD simulators. + if (network->getModules().size() <= 0 ) { + throw std::runtime_error("There are no CFD simulators defined for the Hybrid simulation."); + } - // Initialization of CFD domains - while (! allConverged) { - allConverged = conductCFDSimulation(cfdSimulators, 1); - } + bool allConverged = false; + bool pressureConverged = false; - while (! allConverged || !pressureConverged) { - //std::cout << "######################## Simulation Iteration no. " << iter << " ####################" << std::endl; + // Initialization of CFD domains + while (! allConverged) { + allConverged = conductCFDSimulation(cfdSimulators, 1); + } - // conduct CFD simulations - //std::cout << "[Simulation] Conduct CFD simulation " << iter <<"..." << std::endl; - allConverged = conductCFDSimulation(cfdSimulators, 10); - - // compute nodal analysis again - //std::cout << "[Simulation] Conduct nodal analysis " << iter <<"..." << std::endl; - pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); + while (! allConverged || !pressureConverged) { + // conduct CFD simulations + allConverged = conductCFDSimulation(cfdSimulators, 10); + // compute nodal analysis again + pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); - } + } - #ifdef VERBOSE - if (pressureConverged && allConverged) { - std::cout << "[Simulation] All pressures have converged." << std::endl; - } - printResults(); - #endif + #ifdef VERBOSE + if (pressureConverged && allConverged) { + std::cout << "[Simulation] All pressures have converged." << std::endl; + } + printResults(); + #endif + + saveState(); + } + + // Mixing Hybrid simulation + if (this->simType == Type::Hybrid && this->platform == Platform::Mixing) { + + // Catch runtime error, not enough CFD simulators. + if (network->getModules().size() <= 0 ) { + throw std::runtime_error("There are no CFD simulators defined for the Hybrid simulation."); + } + + bool allConverged = false; + bool pressureConverged = false; + + // Initialization of NS CFD domains + while (! allConverged) { + allConverged = conductCFDSimulation(cfdSimulators, 1); + } + + // Obtain overal steady-state flow result + while (! allConverged || !pressureConverged) { + // conduct CFD simulations + allConverged = conductCFDSimulation(cfdSimulators, 10); + // compute nodal analysis again + pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); + } + + #ifdef VERBOSE + printResults(); + std::cout << "[Simulation] All pressures have converged." << std::endl; + #endif + saveState(); + + // Couple the resulting CFD flow field to the AD fields + coupleNsAdLattices(cfdSimulators); + + // Obtain overal steady-state concentration results + bool concentrationConverged = false; + while (!concentrationConverged) { + concentrationConverged = conductADSimulation(cfdSimulators); + this->mixingModel->propagateSpecies(network, this); + } + } + + // OoC Hybrid simulation + if (this->simType == Type::Hybrid && this->platform == Platform::Ooc) { + + // Catch runtime error, not enough CFD simulators. + if (network->getModules().size() <= 0 ) { + throw std::runtime_error("There are no CFD simulators defined for the Hybrid simulation."); } + + bool allConverged = false; + bool pressureConverged = false; + + // Initialization of CFD domains + while (! allConverged) { + allConverged = conductCFDSimulation(cfdSimulators, 1); + } + + while (! allConverged || !pressureConverged) { + // conduct CFD simulations + allConverged = conductCFDSimulation(cfdSimulators, 10); + // compute nodal analysis again + pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); + } + + #ifdef VERBOSE + if (pressureConverged && allConverged) { + std::cout << "[Simulation] All pressures have converged." << std::endl; + } + printResults(); + #endif + saveState(); } @@ -607,14 +742,15 @@ namespace sim { */ if (simType == Type::Abstract && platform == Platform::Mixing) { T timestep = 0.0; + + // compute nodal analysis + nodalAnalysis->conductNodalAnalysis(); while(true) { if (iteration >= 1000) { throw "Max iterations exceeded."; break; } - // compute nodal analysis - nodalAnalysis->conductNodalAnalysis(); // Update and propagate the mixtures if (this->mixingModel->isInstantaneous()){ @@ -738,6 +874,64 @@ namespace sim { cfdSimulator->prepareLattice(); } } + + if (this->simType == Type::Hybrid && this->platform == Platform::Mixing) { + + #ifdef VERBOSE + std::cout << "[Simulation] Initialize CFD simulators..." << std::endl; + #endif + + // Initialize the CFD simulators + for (auto& [key, cfdSimulator] : cfdSimulators) { + cfdSimulator->lbmInit(fluids[continuousPhase]->getViscosity(), + fluids[continuousPhase]->getDensity()); + } + + // compute nodal analysis + #ifdef VERBOSE + std::cout << "[Simulation] Conduct initial nodal analysis..." << std::endl; + #endif + nodalAnalysis->conductNodalAnalysis(cfdSimulators); + + // Prepare CFD geometry and lattice + #ifdef VERBOSE + std::cout << "[Simulation] Prepare CFD geometry and lattice..." << std::endl; + #endif + + for (auto& [key, cfdSimulator] : cfdSimulators) { + cfdSimulator->prepareGeometry(); + cfdSimulator->prepareLattice(); + } + } + + if (this->simType == Type::Hybrid && this->platform == Platform::Ooc) { + + #ifdef VERBOSE + std::cout << "[Simulation] Initialize CFD simulators..." << std::endl; + #endif + + // Initialize the CFD simulators + for (auto& [key, cfdSimulator] : cfdSimulators) { + cfdSimulator->lbmInit(fluids[continuousPhase]->getViscosity(), + fluids[continuousPhase]->getDensity()); + } + + // compute nodal analysis + #ifdef VERBOSE + std::cout << "[Simulation] Conduct initial nodal analysis..." << std::endl; + #endif + nodalAnalysis->conductNodalAnalysis(cfdSimulators); + + // Prepare CFD geometry and lattice + #ifdef VERBOSE + std::cout << "[Simulation] Prepare CFD geometry and lattice..." << std::endl; + #endif + + for (auto& [key, cfdSimulator] : cfdSimulators) { + cfdSimulator->prepareGeometry(); + cfdSimulator->prepareLattice(); + } + } } template diff --git a/src/simulation/Tissue.h b/src/simulation/Tissue.h new file mode 100644 index 0000000..8c6357b --- /dev/null +++ b/src/simulation/Tissue.h @@ -0,0 +1,90 @@ +/** + * @file Tissue.h + */ + +#pragma once + +#include +#include + +namespace sim { + +// Forward declared dependencies +template +class Specie; + +/** + * @brief Class that describes a tissue. +*/ +template +class Tissue { +private: + + int id; + std::string name; + + std::unordered_map*> species; + std::unordered_map vMax; // vMax values for the enzymatic kinetics per species for this tissue, according to Michaelis-Menten kinetics + std::unordered_map kM; // kM values for the enzymatic kinetics per species for this tissue, according to Michaelis-Menten kinetics + +public: + /** + * @brief Construct a new tissue. + * @param[in] id The id of this tissue. + */ + Tissue(int id); + + /** + * @brief Construct a new tissue. + * @param[in] id The id of this tissue. + * @param[in] specie A pointer to a species that interacts in enzymatic transformation. + * @param[in] vMax The maximum elimination rate of a species in enzymatic transformation. + * @param[in] kM The concentration associated with 0.5vMax for a species in enzymatic transformation. + */ + Tissue(int id, Specie* specie, T vMax, T kM); + + /** + * @brief Construct a new tissue. + * @param[in] id The id of this tissue. + * @param[in] species An unordered map of pointers to species that interact in enzymatic transformation. + * @param[in] vMax An unordered map of the maximum eleimination rate of species in enzymatic transformation. + * @param[in] kM An unordered map of the concentration associated with 0.5vMax for species in enzymatic transformation. + */ + Tissue(int id, std::unordered_map*> species, std::unordered_map vMax, std::unordered_map kM); + + /** + * @brief Get the name of this tissue. + */ + std::string getName(); + + /** + * @brief Set the name of this tissue. + * @param[in] name The new name of this tissue. + */ + void setName(std::string name); + + /** + * @brief Add a new specie interaction to this tissue. + * @param[in] specie A pointer to a species that interacts in enzymatic transformation. + * @param[in] vMax The maximum eleimination rate of a species in enzymatic transformation. + * @param[in] kM The concentration associated with 0.5vMax for a species in enzymatic transformation. + */ + void addSpecie(Specie* species, T vMax, T kM); + + /** + * @brief Add a set of species interactions to this tissue. + * @param[in] species An unordered map of pointers to species that interact in enzymatic transformation. + * @param[in] vMax An unordered map of the maximum eleimination rate of species in enzymatic transformation. + * @param[in] kM An unordered map of the concentration associated with 0.5vMax for species in enzymatic transformation. + */ + void addSpecies(std::unordered_map*> species, std::unordered_map vMax, std::unordered_map kM); + + std::unordered_map*>& getSpecies(); + + T* getVmax(int speciesId); + + T* getkM(int speciesId); + +}; + +} /// namespace sim \ No newline at end of file diff --git a/src/simulation/Tissue.hh b/src/simulation/Tissue.hh new file mode 100644 index 0000000..976ca41 --- /dev/null +++ b/src/simulation/Tissue.hh @@ -0,0 +1,80 @@ +#include "Tissue.h" + +namespace sim { + +template +Tissue::Tissue(int id_) : id(id_) { } + +template +Tissue::Tissue(int id_, Specie* specie_, T vMax_, T kM_) : + id(id_), species(specie_), vMax(vMax_), kM(kM_) { } + +template +Tissue::Tissue(int id_, std::unordered_map*> species_, std::unordered_map vMax_, std::unordered_map kM_) : + id(id_), species(species_), vMax(vMax_), kM(kM_) { } + +template +void Tissue::setName(std::string name_) { + this->name = name_; +} + +template +std::string Tissue::getName() { + return this->name; +} + +template +void Tissue::addSpecie(Specie* specie_, T vMax_, T kM_) { + int specieId = specie_->getId(); + + auto const result = species->try_emplace(specieId, specie_); + if (!result.second) { + std::string errorString = "Cannot add species " + std::to_string(specieId) + " to tissue " + std::to_string(id) + ": This species is already defined for the tissue."; + throw std::invalid_argument(errorString); + } + vMax->try_emplace(specieId, vMax_); + kM->try_emplace(specieId, kM_); +} + +template +void Tissue::addSpecies(std::unordered_map*> species_, std::unordered_map vMax_, std::unordered_map kM_) { + + for (auto& [key, specie_] : species_) { + auto const result = species->try_emplace(key, specie_); + if (!result.second) { + std::string errorString = "Cannot add species " + std::to_string(key) + " to tissue " + std::to_string(id) + ": This species is already defined for the tissue."; + throw std::invalid_argument(errorString); + } + if (vMax_.count(key)) { + vMax->try_emplace(key, vMax_.at(key)); + } else { + std::string errorString = "No vMax value was given for species " + std::to_string(key) + "."; + throw std::invalid_argument(errorString); + } + if (kM_.count(key)) { + kM->try_emplace(key, kM_.at(key)); + } else { + std::string errorString = "No kM value was given for species " + std::to_string(key) + "."; + throw std::invalid_argument(errorString); + } + } +} + +template +std::unordered_map*>& Tissue::getSpecies() { + return species; +} + +template +T* Tissue::getVmax(int speciesId) { + auto vMaxIter = vMax.find(speciesId); + return &(vMaxIter->second); +} + +template +T* Tissue::getkM(int speciesId) { + auto kMIter = kM.find(speciesId); + return &(kMIter->second); +} + +} // namespace sim \ No newline at end of file diff --git a/src/simulation/simulators/CMakeLists.txt b/src/simulation/simulators/CMakeLists.txt index 25d3791..64637a2 100644 --- a/src/simulation/simulators/CMakeLists.txt +++ b/src/simulation/simulators/CMakeLists.txt @@ -1,15 +1,20 @@ set(SOURCE_LIST cfdSimulator.hh olbContinuous.hh + olbMixing.hh + olbOoc.hh ) set(HEADER_LIST cfdSimulator.h olbContinuous.h + olbMixing.h + olbOoc.h ) if(USE_ESSLBM) target_sources(${TARGET_NAME} PUBLIC essContinuous.hh essContinuous.h) + target_sources(${TARGET_NAME} PUBLIC essMixing.hh essMixing.h) endif() target_sources(${TARGET_NAME} PRIVATE ${SOURCE_LIST} ${HEADER_LIST}) diff --git a/src/simulation/simulators/cfdSimulator.h b/src/simulation/simulators/cfdSimulator.h index dd1659d..9c5b885 100755 --- a/src/simulation/simulators/cfdSimulator.h +++ b/src/simulation/simulators/cfdSimulator.h @@ -45,6 +45,11 @@ class CFDSimulator { T alpha; ///< Relaxation factor for convergence between 1D and CFD simulation. + /** + * @brief Define and prepare the coupling of the NS lattice with the AD lattices. + */ + virtual void executeCoupling() { }; + public: CFDSimulator(int id, std::string name, std::string stlFile, std::shared_ptr> cfdModule, std::unordered_map> openings, T alpha, ResistanceModel* ResistanceModel); @@ -119,14 +124,18 @@ class CFDSimulator { */ virtual bool hasConverged() const = 0; - virtual void setPressures(std::unordered_map pressure) = 0; + virtual void storePressures(std::unordered_map pressure) = 0; virtual std::unordered_map getPressures() const = 0; - virtual void setFlowRates(std::unordered_map flowRate) = 0; + virtual void storeFlowRates(std::unordered_map flowRate) = 0; virtual std::unordered_map getFlowRates() const = 0; + virtual void storeConcentrations(std::unordered_map> concentrations) { } + + virtual std::unordered_map> getConcentrations() const { return std::unordered_map>(); } + virtual void setBoundaryValues(int iT) = 0; // virtual functions @@ -135,9 +144,23 @@ class CFDSimulator { virtual void prepareLattice() {} - virtual void writeVTK (int iT) {}; + /** + * @brief Conducts the collide and stream operations of the NS lattice. + */ + virtual void nsSolve() {} + + /** + * @brief Conducts the collide and stream operations of the AD lattice(s). + */ + virtual void adSolve() {} + + virtual void writeVTK (int iT) {} - virtual void getResults (int iT) {}; + virtual void storeCfdResults (int iT) {} + + virtual bool hasAdConverged() const { return false; } + + friend void coupleNsAdLattices(const std::unordered_map>>& cfdSimulators); }; diff --git a/src/simulation/simulators/essContinuous.h b/src/simulation/simulators/essContinuous.h index 4ccce18..eb500b1 100644 --- a/src/simulation/simulators/essContinuous.h +++ b/src/simulation/simulators/essContinuous.h @@ -94,7 +94,7 @@ namespace sim { * @brief Update the values at the module nodes based on the simulation result after stepIter iterations. * @param[in] iT Iteration step. */ - void getResults(); + void storeCfdResults(); /** * @brief Set the pressures at the nodes on the module boundary. diff --git a/src/simulation/simulators/essContinuous.hh b/src/simulation/simulators/essContinuous.hh index 189014d..ad9e29b 100644 --- a/src/simulation/simulators/essContinuous.hh +++ b/src/simulation/simulators/essContinuous.hh @@ -24,7 +24,7 @@ namespace sim{ } template - void essLbmSimulator::getResults() + void essLbmSimulator::storeCfdResults() { for(auto& [key,value] : solver_->getPressures()) pressures[key] = value; @@ -86,7 +86,7 @@ namespace sim{ void essLbmSimulator::solve() { solver_->solve(10, 10, 10); - getResults(); + storeCfdResults(); } template diff --git a/src/simulation/simulators/essMixing.h b/src/simulation/simulators/essMixing.h new file mode 100644 index 0000000..eb500b1 --- /dev/null +++ b/src/simulation/simulators/essMixing.h @@ -0,0 +1,148 @@ +/** + * @file essContinuous.h +*/ + +#pragma once + +#define M_PI 3.14159265358979323846 + +#include +#include +#include +#include +#include + +#include "Module.h" +#include "ModuleOpening.h" +#include "Node.h" +#include "Channel.h" + +#include + +namespace arch { + +// Forward declared dependencies +template +class Module; +template +class Network; +template +class Node; +template +class Opening; + +} + +namespace sim { + + /** + * @brief Class that defines the lbm module which is the interface between the 1D solver and OLB. + */ + template + class essLbmSimulator : public CFDSimulator { + private: + int step = 0; ///< Iteration step of this module. + int stepIter = 1000; ///< Number of iterations for the value tracer. + int maxIter = 1e7; ///< Maximum total iterations. + int theta = 10; ///< Number of OLB iterations per communication iteration. + + bool isConverged = false; ///< Has the module converged? + + T charPhysLength; ///< Characteristic physical length (= width, usually). + T charPhysVelocity; ///< Characteristic physical velocity (expected maximal velocity). + T alpha; ///< Relaxation factor. + T resolution; ///< Resolution of the CFD domain. Gridpoints in charPhysLength. + T epsilon; ///< Convergence criterion. + T relaxationTime; + + std::unordered_map flowRates; ///< Map of fluxes at module nodes. + std::unordered_map pressures; ///< Map of mean pressure values at module nodes. + + std::shared_ptr solver_; + + public: + /** + * @brief Constructor of an lbm module. + * @param[in] id Id of the module. + * @param[in] name Name of the module. + * @param[in] pos Absolute position of the module in _m_, from the bottom left corner of the microfluidic device. + * @param[in] size Size of the module in _m_. + * @param[in] nodes Map of nodes that are on the boundary of the module. + * @param[in] openings Map of the in-/outlets of the module. + */ + essLbmSimulator(int id_, std::string name_, std::string stlFile_, std::shared_ptr> cfdModule, std::unordered_map> openings_, + ResistanceModel* resistanceModel, T charPhysLength_, T charPhysVelocity_, T alpha, T resolution_, T epsilon_, T relaxationTime_); + /** + * @brief Initialize an instance of the LBM solver for this module. + * @param[in] dynViscosity Dynamic viscosity of the simulated fluid in _kg / m s_. + * @param[in] density Density of the simulated fluid in _kg / m^3_. + */ + void lbmInit(T dynViscosity, T density) override; + + /** + * @brief Set the boundary values on the lattice at the module nodes. + * @param[in] iT Iteration step. + */ + void setBoundaryValues(int iT); + + /** + * @brief Conducts the collide and stream operations of the lattice. + */ + void solve() override; + + /** + * @brief Update the values at the module nodes based on the simulation result after stepIter iterations. + * @param[in] iT Iteration step. + */ + void storeCfdResults(); + + /** + * @brief Set the pressures at the nodes on the module boundary. + * @param[in] pressure Map of pressures and node ids. + */ + void setPressures(std::unordered_map pressure) override; + + /** + * @brief Set the flow rates at the nodes on the module boundary. + * @param[in] flowRate Map of flow rates and node ids. + */ + void setFlowRates(std::unordered_map flowRate) override; + + /** + * @brief Get the pressures at the boundary nodes. + * @returns Pressures in Pa. + */ + std::unordered_map getPressures() const override; + + /** + * @brief Get the flow rates at the boundary nodes. + * @returns Flow rates in m^3/s. + */ + std::unordered_map getFlowRates() const override; + + /** + * @brief Get the ground nodes of the module. + * @returns Ground nodes. + */ + std::unordered_map getGroundNodes() { + return this->groundNodes; + } + + /** + * @brief Get the number of iterations for the value tracer. + * @returns Number of iterations for the value tracer. + */ + int getStepIter() const { + return stepIter; + } + + /** + * @brief Returns whether the module has converged or not. + * @returns Boolean for module convergence. + */ + bool hasConverged() const override; + + void setVtkFolder(std::string vtkFolder); + + }; +} diff --git a/src/simulation/simulators/essMixing.hh b/src/simulation/simulators/essMixing.hh new file mode 100644 index 0000000..28c0319 --- /dev/null +++ b/src/simulation/simulators/essMixing.hh @@ -0,0 +1,135 @@ +#include "essContinuous.h" +#include + +#include +#include + +namespace sim{ + + template + essLbmSimulator::essLbmSimulator(int id_, std::string name_, std::string stlFile_, std::shared_ptr> cfdModule_, std::unordered_map> openings_, + ResistanceModel* resistanceModel_, T charPhysLength_, T charPhysVelocity_, T alpha_, T resolution_, T epsilon_, T relaxationTime_) : + CFDSimulator(id_, name_, stlFile_, cfdModule_, openings_, alpha_, resistanceModel_), + charPhysLength(charPhysLength_), charPhysVelocity(charPhysVelocity_), resolution(resolution_), + epsilon(epsilon_), relaxationTime(relaxationTime_) + { + this->cfdModule->setModuleTypeEssLbm(); + allNodes = cfdModule_->getNodes(); + } + + template + void essLbmSimulator::setBoundaryValues(int iT) + { + setFlowRates(flowRates); + setPressures(pressures); + } + + template + void essLbmSimulator::storeCfdResults() + { + for(auto& [key,value] : solver_->getPressures()) + pressures[key] = value; + for(auto& [key,value] : solver_->getFlowRates()) + flowRates[key] = value; + } + + template + void essLbmSimulator::lbmInit(T dynViscosity, T density) + { + + std::string work_dir = "/home/michel/Git/mmft-hybrid-simulator/build/"; + const auto& allNodes = this->moduleNetwork->getNodes(); + std::unordered_map nodes(allNodes.size()); + std::unordered_map openings; + + for(auto& op : this->moduleOpenings) + { + ess::Opening opening; + opening.width = op.second.width; + opening.height = op.second.height; + opening.normal = {op.second.normal[0],op.second.normal[1],op.second.normal[2]}; + openings.try_emplace(op.first, opening); + + ess::BoundaryNode node(allNodes.at(op.first)->getPosition()[0], + allNodes.at(op.first)->getPosition()[1], + allNodes.at(op.first)->getPosition()[2], + allNodes.at(op.first)->getGround()); + + nodes.try_emplace(op.first, node); + } + + solver_ = std::make_shared(work_dir, this->stlFile, openings, nodes, charPhysLength, charPhysVelocity, 1.0f / resolution, epsilon, relaxationTime, density, dynViscosity); + solver_->prepareLattice(); + + #ifdef DEBUG + // There must be more than 1 node to have meaningful flow in the module domain + assert(this->moduleOpenings.size() > 1); + // We must have exactly one opening assigned to each boundaryNode + assert(this->moduleOpenings.size() == this->cfdModule->getNodes().size()); + #endif + + // Initialize pressure, flowRate and resistance value-containers + for (auto& [key, node] : this->moduleOpenings) + { + pressures.try_emplace(key, 0.0f); + flowRates.try_emplace(key, 0.0f); + } + + setFlowRates(flowRates); + setPressures(pressures); + + #ifdef VERBOSE + std::cout << "[essLbmModule] lbmInit " << this->name << "... OK" << std::endl; + #endif + } + + template + void essLbmSimulator::solve() + { + solver_->solve(10, 10, 10); + storeCfdResults(); + } + + template + void essLbmSimulator::setPressures(std::unordered_map pressure_) + { + std::unordered_map interface; + for(auto& [key, value] : pressure_) + interface.try_emplace(key, value); + pressures = pressure_; + solver_->setPressures(interface); + } + + template + void essLbmSimulator::setFlowRates(std::unordered_map flowRate_) + { + std::unordered_map interface; + for(auto& [key, value] : flowRate_) + interface.try_emplace(key, value); + flowRates = flowRate_; + solver_->setFlowRates(interface); + } + + template + std::unordered_map essLbmSimulator::getPressures() const { + return pressures; + } + + + template + std::unordered_map essLbmSimulator::getFlowRates() const { + return flowRates; + } + + template + bool essLbmSimulator::hasConverged() const + { + return solver_->hasConverged(); + } + + template + void essLbmSimulator::setVtkFolder(std::string vtkFolder_) { + this->vtkFolder = vtkFolder_; + } + +} // namespace arch diff --git a/src/simulation/simulators/olbContinuous.h b/src/simulation/simulators/olbContinuous.h index c9e23fb..f35c11c 100644 --- a/src/simulation/simulators/olbContinuous.h +++ b/src/simulation/simulators/olbContinuous.h @@ -44,7 +44,7 @@ using NoDynamics = olb::NoDynamics; using BGKdynamics = olb::BGKdynamics; using BounceBack = olb::BounceBack; -private: +protected: int step = 0; ///< Iteration step of this module. int stepIter = 1000; ///< Number of iterations for the value tracer. int maxIter = 1e7; ///< Maximum total iterations. @@ -87,6 +87,36 @@ using BounceBack = olb::BounceBack; return *lattice; } + void setOutputDir(); + + virtual void initValueContainers(); + + void initNsConverter(T dynViscosity, T density); + + void initNsConvergeTracker(); + + virtual void prepareNsLattice(const T omega); + + void initPressureIntegralPlane(); + + void initFlowRateIntegralPlane(); + + void initNsLattice(const T omega); + + void readGeometryStl(const bool print); + + void readOpenings(); + + void setFlowProfile2D(int key, T openingWidth); + + void setPressure2D(int key); + + /** + * @brief Update the values at the module nodes based on the simulation result after stepIter iterations. + * @param[in] iT Iteration step. + */ + void storeCfdResults(int iT); + public: /** * @brief Constructor of an lbm module. @@ -135,12 +165,6 @@ using BounceBack = olb::BounceBack; */ void solve(); - /** - * @brief Update the values at the module nodes based on the simulation result after stepIter iterations. - * @param[in] iT Iteration step. - */ - void getResults(int iT); - /** * @brief Write the vtk file with results of the CFD simulation to file system. * @param[in] iT Iteration step. @@ -148,16 +172,16 @@ using BounceBack = olb::BounceBack; void writeVTK(int iT); /** - * @brief Set the pressures at the nodes on the module boundary. + * @brief Store the abstract pressures at the nodes on the module boundary in the simulator. * @param[in] pressure Map of pressures and node ids. */ - void setPressures(std::unordered_map pressure); + void storePressures(std::unordered_map pressure); /** - * @brief Set the flow rates at the nodes on the module boundary. + * @brief Store the abstract flow rates at the nodes on the module boundary in the simulator. * @param[in] flowRate Map of flow rates and node ids. */ - void setFlowRates(std::unordered_map flowRate); + void storeFlowRates(std::unordered_map flowRate); /** * @brief Get the pressures at the boundary nodes. diff --git a/src/simulation/simulators/olbContinuous.hh b/src/simulation/simulators/olbContinuous.hh index 0f5ba9d..be728ae 100644 --- a/src/simulation/simulators/olbContinuous.hh +++ b/src/simulation/simulators/olbContinuous.hh @@ -3,16 +3,6 @@ namespace sim{ -template -void lbmSimulator::setPressures(std::unordered_map pressure_) { - this->pressures = pressure_; -} - -template -void lbmSimulator::setFlowRates(std::unordered_map flowRate_) { - this->flowRates = flowRate_; -} - template lbmSimulator::lbmSimulator ( int id_, std::string name_, std::string stlFile_, std::shared_ptr> cfdModule_, std::unordered_map> openings_, @@ -20,84 +10,191 @@ lbmSimulator::lbmSimulator ( CFDSimulator(id_, name_, stlFile_, cfdModule_, openings_, alpha_, resistanceModel_), charPhysLength(charPhysLength_), charPhysVelocity(charPhysVelocity_), resolution(resolution_), epsilon(epsilon_), relaxationTime(relaxationTime_) - { - this->cfdModule->setModuleTypeLbm(); - } +{ + this->cfdModule->setModuleTypeLbm(); +} + +template +void lbmSimulator::lbmInit (T dynViscosity, T density) { + + setOutputDir(); + initValueContainers(); + initNsConverter(dynViscosity, density); + initNsConvergeTracker(); + + #ifdef VERBOSE + std::cout << "[lbmSimulator] lbmInit " << this->name << "... OK" << std::endl; + #endif +} template void lbmSimulator::prepareGeometry () { bool print = false; - stlReader = std::make_shared>(this->stlFile, converter->getConversionFactorLength()); #ifdef VERBOSE print = true; - std::cout << "[lbmModule] reading STL file " << this->name << "... OK" << std::endl; #endif - stl2Dindicator = std::make_shared>(*stlReader); + + readGeometryStl(print); + readOpenings(); + + this->geometry->clean(print); + this->geometry->checkForErrors(print); + + if (print) { + std::cout << "[lbmSimulator] prepare geometry " << this->name << "... OK" << std::endl; + } +} + +template +void lbmSimulator::prepareLattice () { + + const T omega = converter->getLatticeRelaxationFrequency(); + + prepareNsLattice(omega); + initPressureIntegralPlane(); + initFlowRateIntegralPlane(); + initNsLattice(omega); + #ifdef VERBOSE - std::cout << "[lbmModule] create 2D indicator " << this->name << "... OK" << std::endl; + std::cout << "[lbmSimulator] prepare lattice " << this->name << "... OK" << std::endl; #endif +} - olb::Vector origin(-0.5*converter->getConversionFactorLength(), -0.5*converter->getConversionFactorLength()); - olb::Vector extend(this->cfdModule->getSize()[0] + converter->getConversionFactorLength(), this->cfdModule->getSize()[1] + converter->getConversionFactorLength()); - olb::IndicatorCuboid2D cuboid(extend, origin); - cuboidGeometry = std::make_shared> (cuboid, converter->getConversionFactorLength(), 1); - loadBalancer = std::make_shared> (*cuboidGeometry); - geometry = std::make_shared> ( - *cuboidGeometry, - *loadBalancer); +template +void lbmSimulator::setBoundaryValues (int iT) { - #ifdef VERBOSE - std::cout << "[lbmModule] generate geometry " << this->name << "... OK" << std::endl; - #endif + for (auto& [key, Opening] : this->moduleOpenings) { + if (this->groundNodes.at(key)) { + setFlowProfile2D(key, Opening.width); + } else { + setPressure2D(key); + } + } +} - this->geometry->rename(0, 2); - this->geometry->rename(2, 1, *stl2Dindicator); +template +void lbmSimulator::writeVTK (int iT) { + bool print = false; #ifdef VERBOSE - std::cout << "[lbmModule] generate 2D geometry from STL " << this->name << "... OK" << std::endl; + print = true; #endif - for (auto& [key, Opening] : this->moduleOpenings ) { - // The unit vector pointing to the extend (opposite origin) of the opening - T x_origin = Opening.node->getPosition()[0] - this->cfdModule->getPosition()[0] - - 0.5*Opening.width*Opening.tangent[0]; - T y_origin = Opening.node->getPosition()[1] - this->cfdModule->getPosition()[1] - - 0.5*Opening.width*Opening.tangent[1]; + olb::SuperVTMwriter2D vtmWriter( this->name ); + // Writes geometry to file system + if (iT == 0) { + olb::SuperLatticeGeometry2D writeGeometry (getLattice(), getGeometry()); + vtmWriter.write(writeGeometry); + vtmWriter.createMasterFile(); + this->vtkFile = olb::singleton::directories().getVtkOutDir() + olb::createFileName( this->name ) + ".pvd"; + } + + if (iT % 1000 == 0) { - // The unit vector pointing to the extend - T x_extend = Opening.width*Opening.tangent[0] - converter->getConversionFactorLength()*Opening.normal[0]; - T y_extend = Opening.width*Opening.tangent[1] - converter->getConversionFactorLength()*Opening.normal[1]; + olb::SuperLatticePhysVelocity2D velocity(getLattice(), getConverter()); + olb::SuperLatticePhysPressure2D pressure(getLattice(), getConverter()); + olb::SuperLatticeDensity2D latDensity(getLattice()); + vtmWriter.addFunctor(velocity); + vtmWriter.addFunctor(pressure); + vtmWriter.addFunctor(latDensity); + + // write vtk to file system + vtmWriter.write(iT); + this->vtkFile = olb::singleton::directories().getVtkOutDir() + "data/" + olb::createFileName( this->name, iT ) + ".vtm"; + converge->takeValue(getLattice().getStatistics().getAverageEnergy(), print); + } + if (iT %1000 == 0) { + #ifdef VERBOSE + std::cout << "[writeVTK] " << this->name << " currently at timestep " << iT << std::endl; + for (auto& [key, Opening] : this->moduleOpenings) { + if (this->groundNodes.at(key)) { + meanPressures.at(key)->print(); + } else { + fluxes.at(key)->print(); + } + } + #endif + } - // Extend can only have positive values, hence the following transformation - if (x_extend < 0 ){ - x_extend *= -1; - x_origin -= x_extend; - } - if (y_extend < 0 ){ - y_extend *= -1; - y_origin -= y_extend; + converge->takeValue(getLattice().getStatistics().getAverageEnergy(), print); + + if (iT%100 == 0) { + if (converge->hasConverged()) { + isConverged = true; } + } - olb::Vector originO (x_origin, y_origin); - olb::Vector extendO (x_extend, y_extend); - olb::IndicatorCuboid2D opening(extendO, originO); - - this->geometry->rename(2, key+3, 1, opening); +} + +template +void lbmSimulator::solve() { + // theta = 10 + this->setBoundaryValues(step); + for (int iT = 0; iT < 10; ++iT){ + writeVTK(step); + lattice->collideAndStream(); + step += 1; } + storeCfdResults(step); +} - this->geometry->clean(print); - this->geometry->checkForErrors(print); +template +void lbmSimulator::setOutputDir () { + if (!std::filesystem::is_directory(this->vtkFolder) || !std::filesystem::exists(this->vtkFolder)) { + std::filesystem::create_directory(this->vtkFolder); + } + + olb::singleton::directories().setOutputDir( this->vtkFolder+"/" ); // set output directory +} + +template +void lbmSimulator::initValueContainers () { + + // Initialize pressure and flow rate value-containers + for (auto& [key, node] : this->moduleOpenings) { + pressures.try_emplace(key, (T) 0.0); + flowRates.try_emplace(key, (T) 0.0); + } + +} + +template +void lbmSimulator::initNsConverter (T dynViscosity, T density) { + + T kinViscosity = dynViscosity/density; + + #ifdef DEBUG + // There must be more than 1 node to have meaningful flow in the module domain + assert(this->moduleOpenings.size() > 1); + // We must have exactly one opening assigned to each boundaryNode + assert(this->moduleOpenings.size() == this->cfdModule->getNodes().size()); + #endif + + this->converter = std::make_shared>( + resolution, + relaxationTime, + charPhysLength, + charPhysVelocity, + kinViscosity, + density + ); #ifdef VERBOSE - std::cout << "[lbmModule] prepare geometry " << this->name << "... OK" << std::endl; + this->converter->print(); #endif + } template -void lbmSimulator::prepareLattice () { - const T omega = converter->getLatticeRelaxationFrequency(); +void lbmSimulator::initNsConvergeTracker () { + // Initialize a convergence tracker for pressure + this->converge = std::make_unique> (stepIter, epsilon); +} + +template +void lbmSimulator::prepareNsLattice (const T omega) { lattice = std::make_shared>(getGeometry()); @@ -124,6 +221,11 @@ void lbmSimulator::prepareLattice () { } } +} + +template +void lbmSimulator::initPressureIntegralPlane() { + // Initialize the integral fluxes for the in- and outlets for (auto& [key, Opening] : this->moduleOpenings) { @@ -139,7 +241,23 @@ void lbmSimulator::prepareLattice () { position, Opening.tangent, materials); this->meanPressures.try_emplace(key, meanPressure); flowProfiles.try_emplace(key, std::make_shared>(getGeometry(), 0, (T) 0.0, (T) 0.0)); - } else { + } + } +} + +template +void lbmSimulator::initFlowRateIntegralPlane() { + + // Initialize the integral fluxes for the in- and outlets + for (auto& [key, Opening] : this->moduleOpenings) { + + T posX = Opening.node->getPosition()[0] - this->cfdModule->getPosition()[0]; + T posY = Opening.node->getPosition()[1] - this->cfdModule->getPosition()[1]; + + std::vector position = {posX, posY}; + std::vector materials = {1, key+3}; + + if (!this->groundNodes.at(key)) { std::shared_ptr> flux; flux = std::make_shared< olb::SuperPlaneIntegralFluxVelocity2D > (getLattice(), getConverter(), getGeometry(), position, Opening.tangent, materials); @@ -148,174 +266,122 @@ void lbmSimulator::prepareLattice () { } } +} + +template +void lbmSimulator::initNsLattice (const T omega) { // Initialize lattice with relaxation frequency omega lattice->template setParameter(omega); lattice->initialize(); - - #ifdef VERBOSE - std::cout << "[lbmModule] prepare lattice " << this->name << "... OK" << std::endl; - #endif } + template -void lbmSimulator::setBoundaryValues (int iT) { +void lbmSimulator::readGeometryStl (const bool print) { - T pressureLow = -1.0; - for (auto& [key, pressure] : pressures) { - if ( pressureLow < 0.0 ) { - pressureLow = pressure; - } - if ( pressure < pressureLow ) { - pressureLow = pressure; - } - } + stlReader = std::make_shared>(this->stlFile, converter->getConversionFactorLength()); - for (auto& [key, Opening] : this->moduleOpenings) { - if (this->groundNodes.at(key)) { - T maxVelocity = (3./2.)*(flowRates[key]/(Opening.width)); - T distance2Wall = getConverter().getConversionFactorLength()/2.; - this->flowProfiles.at(key) = std::make_shared>(getGeometry(), key+3, getConverter().getLatticeVelocity(maxVelocity), distance2Wall); - getLattice().defineU(getGeometry(), key+3, *this->flowProfiles.at(key)); - } else { - T rhoV = getConverter().getLatticeDensityFromPhysPressure((pressures[key])); - this->densities.at(key) = std::make_shared>(rhoV); - getLattice().defineRho(getGeometry(), key+3, *this->densities.at(key)); - } + if (print) { + std::cout << "[lbmSimulator] reading STL file " << this->name << "... OK" << std::endl; } + + stl2Dindicator = std::make_shared>(*stlReader); -} - -template -void lbmSimulator::getResults (int iT) { - int input[1] = { }; - T output[10]; - - for (auto& [key, Opening] : this->moduleOpenings) { - if (this->groundNodes.at(key)) { - meanPressures.at(key)->operator()(output, input); - T newPressure = output[0]/output[1]; - pressures.at(key) = newPressure; - if (iT % 1000 == 0) { - #ifdef VERBOSE - meanPressures.at(key)->print(); - #endif - } - } else { - fluxes.at(key)->operator()(output,input); - flowRates.at(key) = output[0]; - if (iT % 1000 == 0) { - #ifdef VERBOSE - fluxes.at(key)->print(); - #endif - } - } + if (print) { + std::cout << "[lbmSimulator] create 2D indicator " << this->name << "... OK" << std::endl; } -} -template -void lbmSimulator::lbmInit (T dynViscosity, - T density) { + olb::Vector origin(-0.5*converter->getConversionFactorLength(), -0.5*converter->getConversionFactorLength()); + olb::Vector extend(this->cfdModule->getSize()[0] + converter->getConversionFactorLength(), this->cfdModule->getSize()[1] + converter->getConversionFactorLength()); + olb::IndicatorCuboid2D cuboid(extend, origin); + cuboidGeometry = std::make_shared> (cuboid, converter->getConversionFactorLength(), 1); + loadBalancer = std::make_shared> (*cuboidGeometry); + geometry = std::make_shared> (*cuboidGeometry, *loadBalancer); - if (!std::filesystem::is_directory(this->vtkFolder) || !std::filesystem::exists(this->vtkFolder)) { - std::filesystem::create_directory(this->vtkFolder); + if (print) { + std::cout << "[lbmSimulator] generate geometry " << this->name << "... OK" << std::endl; } - olb::singleton::directories().setOutputDir( this->vtkFolder+"/" ); // set output directory - - T kinViscosity = dynViscosity/density; - - #ifdef DEBUG - // There must be more than 1 node to have meaningful flow in the module domain - assert(this->moduleOpenings.size() > 1); - // We must have exactly one opening assigned to each boundaryNode - assert(this->moduleOpenings.size() == this->cfdModule->getNodes().size()); - #endif - - this->converter = std::make_shared>( - resolution, - relaxationTime, - charPhysLength, - charPhysVelocity, - kinViscosity, - density - ); - - #ifdef VERBOSE - this->converter->print(); - #endif + this->geometry->rename(0, 2); + this->geometry->rename(2, 1, *stl2Dindicator); - // Initialize pressure, flowRate and resistance value-containers - for (auto& [key, node] : this->moduleOpenings) { - pressures.try_emplace(key, (T) 0.0); - flowRates.try_emplace(key, (T) 0.0); + if (print) { + std::cout << "[lbmSimulator] generate 2D geometry from STL " << this->name << "... OK" << std::endl; } - - // Initialize a convergence tracker for pressure - this->converge = std::make_unique> (stepIter, epsilon); - - #ifdef VERBOSE - std::cout << "[lbmModule] lbmInit " << this->name << "... OK" << std::endl; - #endif } template -void lbmSimulator::writeVTK (int iT) { +void lbmSimulator::readOpenings () { - bool print = false; - #ifdef VERBOSE - print = true; - #endif + for (auto& [key, Opening] : this->moduleOpenings ) { + // The unit vector pointing to the extend (opposite origin) of the opening + T x_origin = Opening.node->getPosition()[0] - this->cfdModule->getPosition()[0] + - 0.5*Opening.width*Opening.tangent[0]; + T y_origin = Opening.node->getPosition()[1] - this->cfdModule->getPosition()[1] + - 0.5*Opening.width*Opening.tangent[1]; + + // The unit vector pointing to the extend + T x_extend = Opening.width*Opening.tangent[0] - converter->getConversionFactorLength()*Opening.normal[0]; + T y_extend = Opening.width*Opening.tangent[1] - converter->getConversionFactorLength()*Opening.normal[1]; - olb::SuperVTMwriter2D vtmWriter( this->name ); - // Writes geometry to file system - if (iT == 0) { - olb::SuperLatticeGeometry2D writeGeometry (getLattice(), getGeometry()); - vtmWriter.write(writeGeometry); - vtmWriter.createMasterFile(); - this->vtkFile = olb::singleton::directories().getVtkOutDir() + olb::createFileName( this->name ) + ".pvd"; - } + // Extend can only have positive values, hence the following transformation + if (x_extend < 0 ){ + x_extend *= -1; + x_origin -= x_extend; + } + if (y_extend < 0 ){ + y_extend *= -1; + y_origin -= y_extend; + } - if (iT % 1000 == 0) { - - olb::SuperLatticePhysVelocity2D velocity(getLattice(), getConverter()); - olb::SuperLatticePhysPressure2D pressure(getLattice(), getConverter()); - olb::SuperLatticeDensity2D latDensity(getLattice()); - vtmWriter.addFunctor(velocity); - vtmWriter.addFunctor(pressure); - vtmWriter.addFunctor(latDensity); + olb::Vector originO (x_origin, y_origin); + olb::Vector extendO (x_extend, y_extend); + olb::IndicatorCuboid2D opening(extendO, originO); - // write vtk to file system - vtmWriter.write(iT); - this->vtkFile = olb::singleton::directories().getVtkOutDir() + "data/" + olb::createFileName( this->name, iT ) + ".vtm"; - converge->takeValue(getLattice().getStatistics().getAverageEnergy(), print); - } - if (iT %1000 == 0) { - #ifdef VERBOSE - std::cout << "[writeVTK] " << this->name << " currently at timestep " << iT << std::endl; - #endif + this->geometry->rename(2, key+3, 1, opening); } +} - converge->takeValue(getLattice().getStatistics().getAverageEnergy(), print); +template +void lbmSimulator::storePressures(std::unordered_map pressure_) { + this->pressures = pressure_; +} - if (iT%100 == 0) { - if (converge->hasConverged()) { - isConverged = true; - } - } +template +void lbmSimulator::storeFlowRates(std::unordered_map flowRate_) { + this->flowRates = flowRate_; +} +template +void lbmSimulator::setFlowProfile2D (int openingKey, T openingWidth) { + T maxVelocity = (3./2.)*(flowRates[openingKey]/(openingWidth)); + T distance2Wall = getConverter().getConversionFactorLength()/2.; + this->flowProfiles.at(openingKey) = std::make_shared>(getGeometry(), openingKey+3, getConverter().getLatticeVelocity(maxVelocity), distance2Wall); + getLattice().defineU(getGeometry(), openingKey+3, *this->flowProfiles.at(openingKey)); } template -void lbmSimulator::solve() { - // theta = 10 - for (int iT = 0; iT < 10; ++iT){ - this->setBoundaryValues(step); - writeVTK(step); - lattice->collideAndStream(); - step += 1; - } - getResults(step); +void lbmSimulator::setPressure2D (int openingKey) { + T rhoV = getConverter().getLatticeDensityFromPhysPressure((pressures[openingKey])); + this->densities.at(openingKey) = std::make_shared>(rhoV); + getLattice().defineRho(getGeometry(), openingKey+3, *this->densities.at(openingKey)); } +template +void lbmSimulator::storeCfdResults (int iT) { + int input[1] = { }; + T output[10]; + + for (auto& [key, Opening] : this->moduleOpenings) { + if (this->groundNodes.at(key)) { + meanPressures.at(key)->operator()(output, input); + T newPressure = output[0]/output[1]; + pressures.at(key) = newPressure; + } else { + fluxes.at(key)->operator()(output,input); + flowRates.at(key) = output[0]; + } + } +} } // namespace arch \ No newline at end of file diff --git a/src/simulation/simulators/olbMixing.h b/src/simulation/simulators/olbMixing.h new file mode 100644 index 0000000..77e2d10 --- /dev/null +++ b/src/simulation/simulators/olbMixing.h @@ -0,0 +1,184 @@ +/** + * @file olbMixing.h + */ + +#pragma once + +#define M_PI 3.14159265358979323846 + +#include +#include +#include +#include +#include + +namespace arch { + +// Forward declared dependencies +template +class Module; +template +class Network; +template +class Node; +template +class Opening; + +} + +namespace sim { + +/** + * @brief Class that defines the lbm module which is the interface between the 1D solver and OLB. +*/ +template +class lbmMixingSimulator : public lbmSimulator { + +using DESCRIPTOR = olb::descriptors::D2Q9<>; +using NoDynamics = olb::NoDynamics; +using BGKdynamics = olb::BGKdynamics; +using BounceBack = olb::BounceBack; + +using ADDESCRIPTOR = olb::descriptors::D2Q5; +using ADDynamics = olb::AdvectionDiffusionBGKdynamics; +using NoADDynamics = olb::NoDynamics; + +protected: + std::unordered_map> concentrations; ///< Vector of concentration values at module nodes. > + + std::unordered_map*> species; + + T adRelaxationTime; ///< Relaxation time (tau) for the OLB solver. + + std::unordered_map averageDensities; + std::unordered_map custConverges; + + std::unordered_map>> adLattices; ///< The LBM lattice on the geometry. + std::unordered_map>> adConverges; ///< Value tracer to track convergence. + std::unordered_map>> adConverters; ///< Object that stores conversion factors from phyical to lattice parameters. + + std::unordered_map fluxWall; + T zeroFlux = 0.0; + + std::unordered_map>>> concentrationProfiles; + std::unordered_map>>> meanConcentrations; ///< Map of mean pressure values at module nodes. + + auto& getAdConverter(int key) { + return *adConverters.at(key); + } + + auto& getAdLattice(int key) { + return *adLattices.at(key); + } + + void initValueContainers() override; + + void initAdConverters(T density); + + void initAdConvergenceTracker(); + + void prepareAdLattice(const T omega, int speciesId); + + void initConcentrationIntegralPlane(int adKey); + + void initAdLattice(int adKey); + + /** + * @brief Define and prepare the coupling of the NS lattice with the AD lattices. + */ + void prepareCoupling(); + + /** + * @brief Execute the coupling placed onto the NS lattice. + */ + void executeCoupling() override; + + void setConcentration2D(int key); + + /** + * @brief Update the values at the module nodes based on the simulation result after stepIter iterations. + * @param[in] iT Iteration step. + */ + void storeCfdResults(int iT); + +public: + /** + * @brief Constructor of an lbm module. + * @param[in] id Id of the module. + * @param[in] name Name of the module. + * @param[in] pos Absolute position of the module in _m_, from the bottom left corner of the microfluidic device. + * @param[in] size Size of the module in _m_. + * @param[in] nodes Map of nodes that are on the boundary of the module. + * @param[in] openings Map of the in-/outlets of the module. + * @param[in] stlFile STL file that describes the geometry of the CFD domain. + * @param[in] charPhysLength Characteristic physical length of the geometry of the module in _m_. + * @param[in] charPhysVelocity Characteristic physical velocity of the flow in the module in _m/s_. + * @param[in] alpha Relaxation factor for the iterative updates between the 1D and CFD solvers. + * @param[in] resolution Resolution of the CFD mesh in gridpoints per charPhysLength. + * @param[in] epsilon Convergence criterion for the pressure values at nodes on the boundary of the module. + * @param[in] relaxationTime Relaxation time tau for the LBM solver. + */ + lbmMixingSimulator(int id, std::string name, std::string stlFile, std::shared_ptr> cfdModule, std::unordered_map*> species, + std::unordered_map> openings, ResistanceModel* resistanceModel, T charPhysLenth, T charPhysVelocity, + T alpha, T resolution, T epsilon, T relaxationTime=0.932, T adRelaxationTime=0.932); + + /** + * @brief Initialize an instance of the LBM solver for this module. + * @param[in] dynViscosity Dynamic viscosity of the simulated fluid in _kg / m s_. + * @param[in] density Density of the simulated fluid in _kg / m^3_. + */ + void lbmInit(T dynViscosity, T density) override; + + /** + * @brief Prepare the LBM lattice on the LBM geometry. + */ + void prepareLattice() override; + + /** + * @brief Set the boundary values on the lattice at the module nodes. + * @param[in] iT Iteration step. + */ + void setBoundaryValues(int iT) override; + + /** + * @brief Conducts the collide and stream operations of the lattice. + */ + void solve() override; + + /** + * @brief Conducts the collide and stream operations of the NS lattice. + */ + void nsSolve() override; + + /** + * @brief Conducts the collide and stream operations of the AD lattice(s). + */ + void adSolve() override; + + /** + * @brief Write the vtk file with results of the CFD simulation to file system. + * @param[in] iT Iteration step. + */ + void writeVTK(int iT) override; + + /** + * @brief Store the abstract concentrations at the nodes on the module boundary in the simulator. + * @param[in] concentrations Map of concentrations and node ids. + */ + void storeConcentrations(std::unordered_map> concentrations) override; + + /** + * @brief Get the concentrations at the boundary nodes. + * @returns Concentrations + */ + std::unordered_map> getConcentrations() const override; + + /** + * @brief Returns whether the module has converged or not. + * @returns Boolean for module convergence. + */ + bool hasAdConverged() const override; + +}; + +} // namespace arch diff --git a/src/simulation/simulators/olbMixing.hh b/src/simulation/simulators/olbMixing.hh new file mode 100644 index 0000000..86f33f5 --- /dev/null +++ b/src/simulation/simulators/olbMixing.hh @@ -0,0 +1,403 @@ +#include "olbMixing.h" +#include + +namespace sim{ + +template +lbmMixingSimulator::lbmMixingSimulator ( + int id_, std::string name_, std::string stlFile_, std::shared_ptr> cfdModule_, std::unordered_map*> species_, + std::unordered_map> openings_, ResistanceModel* resistanceModel_, T charPhysLength_, T charPhysVelocity_, + T alpha_, T resolution_, T epsilon_, T relaxationTime_, T adRelaxationTime_) : + lbmSimulator(id_, name_, stlFile_, cfdModule_, openings_, resistanceModel_, charPhysLength_, charPhysVelocity_, alpha_, resolution_, epsilon_, relaxationTime_), + species(species_), adRelaxationTime(adRelaxationTime_) + { + std::cout << "Creating module and setting its type to lbm" << std::endl; + this->cfdModule->setModuleTypeLbm(); + fluxWall.try_emplace(int(0), &zeroFlux); + } + +template +void lbmMixingSimulator::lbmInit (T dynViscosity, T density) { + + this->setOutputDir(); + initValueContainers(); + this->initNsConverter(dynViscosity, density); + initAdConverters(density); + this->initNsConvergeTracker(); + initAdConvergenceTracker(); + + #ifdef VERBOSE + std::cout << "[lbmSimulator] lbmInit " << this->name << "... OK" << std::endl; + #endif +} + +template +void lbmMixingSimulator::prepareLattice () { + + /** + * Prepare the NS lattice + */ + const T omega = this->converter->getLatticeRelaxationFrequency(); + this->prepareNsLattice(omega); + + /** + * Prepare the AD lattices + */ + for (auto& [speciesId, converter] : adConverges) { + const T adOmega = getAdConverter(speciesId).getLatticeRelaxationFrequency(); + prepareAdLattice(adOmega, speciesId); + } + + /** + * Initialize the integral fluxes for the in- and outlets + */ + this->initPressureIntegralPlane(); + this->initFlowRateIntegralPlane(); + for (auto& [adKey, LatticeAD] : adLattices) { + initConcentrationIntegralPlane(adKey); + } + + /** + * Initialize the lattices + */ + this->initNsLattice(omega); + for (auto& [speciesId, converter] : adConverges) { + initAdLattice(speciesId); + } + + std::cout << "[lbmSimulator] prepare lattice " << this->name << "... OK" << std::endl; + + prepareCoupling(); + + std::cout << "[lbmSimulator] prepare coupling " << this->name << "... OK" << std::endl; +} + +template +void lbmMixingSimulator::setBoundaryValues (int iT) { + + for (auto& [key, Opening] : this->moduleOpenings) { + if (this->groundNodes.at(key)) { + this->setFlowProfile2D(key, Opening.width); + } else { + this->setPressure2D(key); + } + setConcentration2D(key); + } +} + +template +void lbmMixingSimulator::storeCfdResults (int iT) { + int input[1] = { }; + T output[10]; + + for (auto& [key, Opening] : this->moduleOpenings) { + if (this->groundNodes.at(key)) { + this->meanPressures.at(key)->operator()(output, input); + T newPressure = output[0]/output[1]; + this->pressures.at(key) = newPressure; + } else { + this->fluxes.at(key)->operator()(output,input); + this->flowRates.at(key) = output[0]; + } + } + + for (auto& [key, Opening] : this->moduleOpenings) { + // If the node is an outflow, write the concentration value + if (this->flowRates.at(key) < 0.0) { + for (auto& [speciesId, adLattice] : adLattices) { + this->meanConcentrations.at(key).at(speciesId)->operator()(output,input); + this->concentrations.at(key).at(speciesId) = output[0]; + if (iT % 1000 == 0) { + this->meanConcentrations.at(key).at(speciesId)->print(); + } + } + } + } + +} + +template +void lbmMixingSimulator::writeVTK (int iT) { + + bool print = false; + #ifdef VERBOSE + print = true; + #endif + + olb::SuperVTMwriter2D vtmWriter( this->name ); + // Writes geometry to file system + if (iT == 0) { + olb::SuperLatticeGeometry2D writeGeometry (this->getLattice(), this->getGeometry()); + vtmWriter.write(writeGeometry); + vtmWriter.createMasterFile(); + this->vtkFile = olb::singleton::directories().getVtkOutDir() + olb::createFileName( this->name ) + ".pvd"; + } + + if (iT % 1000 == 0) { + + olb::SuperLatticePhysVelocity2D velocity(this->getLattice(), this->getConverter()); + olb::SuperLatticePhysPressure2D pressure(this->getLattice(), this->getConverter()); + olb::SuperLatticeDensity2D latDensity(this->getLattice()); + vtmWriter.addFunctor(velocity); + vtmWriter.addFunctor(pressure); + vtmWriter.addFunctor(latDensity); + + vtmWriter.write(iT); + + // write all concentrations + for (auto& [speciesId, adLattice] : adLattices) { + olb::SuperLatticeDensity2D concentration( getAdLattice(speciesId) ); + concentration.getName() = "concentration " + std::to_string(speciesId); + vtmWriter.write(concentration, iT); + } + + // write vtk to file system + this->vtkFile = olb::singleton::directories().getVtkOutDir() + "data/" + olb::createFileName( this->name, iT ) + ".vtm"; + this->converge->takeValue(this->getLattice().getStatistics().getAverageEnergy(), !print); + for (auto& [key, adConverge] : adConverges) { + //adConverge->takeValue(getAdLattice(key).getStatistics().getAverageRho(), print); + T newRho = getAdLattice(key).getStatistics().getAverageRho(); + if (std::abs(averageDensities.at(key) - newRho) < 1e-5) { + custConverges.at(key) = true; + } + averageDensities.at(key) = newRho; + } + } + if (iT %1000 == 0) { + #ifdef VERBOSE + std::cout << "[writeVTK] " << this->name << " currently at timestep " << iT << std::endl; + #endif + } + + this->converge->takeValue(this->getLattice().getStatistics().getAverageEnergy(), print); + + if (iT%100 == 0) { + if (this->converge->hasConverged()) { + this->isConverged = true; + } + } + +} + +template +void lbmMixingSimulator::solve() { + // theta = 10 + this->setBoundaryValues(this->step); + for (int iT = 0; iT < 10; ++iT){ + this->lattice->collideAndStream(); + this->lattice->executeCoupling(); + for (auto& [speciesId, adLattice] : adLattices) { + adLattice->collideAndStream(); + } + writeVTK(this->step); + this->step += 1; + } + storeCfdResults(this->step); +} + +template +void lbmMixingSimulator::executeCoupling() { + this->lattice->executeCoupling(); + std::cout << "[lbmSimulator] Execute NS-AD coupling " << this->name << "... OK" << std::endl; +} + + +template +void lbmMixingSimulator::nsSolve() { + // theta = 10 + this->setBoundaryValues(this->step); + for (int iT = 0; iT < 10; ++iT){ + this->lattice->collideAndStream(); + writeVTK(this->step); + this->step += 1; + } + storeCfdResults(this->step); +} + +template +void lbmMixingSimulator::adSolve() { + // theta = 10 + this->setBoundaryValues(this->step); + for (int iT = 0; iT < 100; ++iT){ + for (auto& [speciesId, adLattice] : adLattices) { + adLattice->collideAndStream(); + } + writeVTK(this->step); + this->step += 1; + } + storeCfdResults(this->step); +} + +template +void lbmMixingSimulator::initValueContainers () { + // Initialize pressure, flowRate and concentration value-containers + for (auto& [key, node] : this->moduleOpenings) { + this->pressures.try_emplace(key, (T) 0.0); + this->flowRates.try_emplace(key, (T) 0.0); + std::unordered_map tmpConcentrations; + for (auto& [speciesId, speciesPtr] : species) { + tmpConcentrations.try_emplace(speciesId, 0.0); + } + this->concentrations.try_emplace(key, tmpConcentrations); + } +} + +template +void lbmMixingSimulator::initAdConverters (T density) { + for (auto& [speciesId, specie] : species) { + std::shared_ptr> tempAD = std::make_shared> ( + this->resolution, + this->relaxationTime, + this->charPhysLength, + this->charPhysVelocity, + specie->getDiffusivity(), + density + ); + #ifdef VERBOSE + tempAD->print(); + #endif + + this->adConverters.try_emplace(speciesId, tempAD); + } +} + +template +void lbmMixingSimulator::initAdConvergenceTracker () { + // Initialize a convergence tracker for concentrations + for (auto& [speciesId, specie] : species) { + this->adConverges.try_emplace(speciesId, std::make_unique> (1000, 1e-1)); + this->averageDensities.try_emplace(speciesId, T(0.0)); + this->custConverges.try_emplace(speciesId, false); + } +} + +template +void lbmMixingSimulator::prepareAdLattice (const T adOmega, int speciesId) { + + std::shared_ptr> adLattice = std::make_shared>(this->getGeometry()); + + // Initial conditions + std::vector zeroVelocity(T(0), T(0)); + olb::AnalyticalConst2D rhoF(1); + olb::AnalyticalConst2D rhoZero(0); + olb::AnalyticalConst2D uZero(zeroVelocity); + + // Set AD lattice dynamics + adLattice->template defineDynamics(this->getGeometry(), 0); + adLattice->template defineDynamics(this->getGeometry(), 1); + adLattice->template defineDynamics(this->getGeometry(), 2); + + // Set initial conditions + adLattice->defineRhoU(this->getGeometry(), 0, rhoZero, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 0, rhoZero, uZero); + adLattice->defineRhoU(this->getGeometry(), 1, rhoF, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 1, rhoF, uZero); + adLattice->defineRhoU(this->getGeometry(), 2, rhoF, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 2, rhoF, uZero); + + adLattice->template defineField(this->getGeometry(), 0, uZero); + adLattice->template defineField(this->getGeometry(), 1, uZero); + adLattice->template defineField(this->getGeometry(), 2, uZero); + + for (auto& [key, Opening] : this->moduleOpenings) { + adLattice->template defineDynamics(this->getGeometry(), key+3); + adLattice->iniEquilibrium(this->getGeometry(), key+3, rhoF, uZero); + adLattice->template defineField(this->getGeometry(), key+3, uZero); + } + + // Add wall boundary + olb::setFunctionalRegularizedHeatFluxBoundary(*adLattice, adOmega, this->getGeometry(), 2, fluxWall.at(0), fluxWall.at(0)); + + adLattices.try_emplace(speciesId, adLattice); + +} + +template +void lbmMixingSimulator::initAdLattice(int adKey) { + const T adOmega = getAdConverter(adKey).getLatticeRelaxationFrequency(); + getAdLattice(adKey).template setParameter(adOmega); + getAdLattice(adKey).template setParameter(1./12.); + getAdLattice(adKey).initialize(); +} + +template +void lbmMixingSimulator::initConcentrationIntegralPlane(int adKey) { + + // Initialize the integral fluxes for the in- and outlets + for (auto& [key, Opening] : this->moduleOpenings) { + + T posX = Opening.node->getPosition()[0] - this->cfdModule->getPosition()[0]; + T posY = Opening.node->getPosition()[1] - this->cfdModule->getPosition()[1]; + + std::vector position = {posX, posY}; + std::vector materials = {1, key+3}; + + std::unordered_map>> meanConcentration; + std::shared_ptr> concentration; + concentration = std::make_shared> (getAdLattice(adKey), getAdConverter(adKey), this->getGeometry(), + position, Opening.tangent, materials); + meanConcentration.try_emplace(adKey, concentration); + this->meanConcentrations.try_emplace(key, meanConcentration); + } +} + +template +void lbmMixingSimulator::prepareCoupling() { + + std::vector*> adLatticesVec; + std::vector velFactors; + for (auto& [speciesId, adLattice] : adLattices) { + adLatticesVec.emplace_back(adLattices[speciesId].get()); + //velFactors.emplace_back(this->converter->getConversionFactorVelocity() / this->adConverters[speciesId]->getConversionFactorVelocity()); + velFactors.emplace_back(1.0); + } + olb::NavierStokesAdvectionDiffusionSingleCouplingGenerator2D coupling(0, this->converter->getLatticeLength(this->cfdModule->getSize()[0]), 0, this->converter->getLatticeLength(this->cfdModule->getSize()[1]), velFactors); + this->lattice->addLatticeCoupling(coupling, adLatticesVec); +} + +template +void lbmMixingSimulator::setConcentration2D (int key) { + // Set the boundary concentrations for inflows and outflows + // If the boundary is an inflow + if (this->flowRates.at(key) >= 0.0) { + for (auto& [speciesId, adLattice] : adLattices) { + olb::setAdvectionDiffusionTemperatureBoundary(*adLattice, getAdConverter(speciesId).getLatticeRelaxationFrequency(), this->getGeometry(), key+3); + olb::AnalyticalConst2D one( 1. ); + olb::AnalyticalConst2D inConc(concentrations.at(key).at(speciesId)); + adLattice->defineRho(this->getGeometry(), key+3, one + inConc); + } + } + // If the boundary is an outflow + else if (this->flowRates.at(key) < 0.0) { + for (auto& [speciesId, adLattice] : adLattices) { + olb::setRegularizedHeatFluxBoundary(*adLattice, getAdConverter(speciesId).getLatticeRelaxationFrequency(), this->getGeometry(), key+3, &zeroFlux); + } + } + else { + std::cerr << "Error: Invalid Flow Type Boundary Node." << std::endl; + exit(1); + } +} + +template +void lbmMixingSimulator::storeConcentrations(std::unordered_map> concentrations_) { + this->concentrations = concentrations_; +} + +template +std::unordered_map> lbmMixingSimulator::getConcentrations() const { + return this->concentrations; +} + +template +bool lbmMixingSimulator::hasAdConverged() const { + bool c = true; + for (auto& [key, converge] : custConverges) { + if (!converge) { + c = false; + } + } + return c; +}; + +} // namespace arch \ No newline at end of file diff --git a/src/simulation/simulators/olbOoc.h b/src/simulation/simulators/olbOoc.h new file mode 100644 index 0000000..6c2b954 --- /dev/null +++ b/src/simulation/simulators/olbOoc.h @@ -0,0 +1,112 @@ +/** + * @file olbOoc.h + */ + +#pragma once + +#define M_PI 3.14159265358979323846 + +#include +#include +#include +#include +#include + +namespace arch { + +// Forward declared dependencies +template +class Module; +template +class Network; +template +class Node; +template +class Opening; + +} + +namespace sim { + +// Forward declared dependencies +template +class Tissue; + +/** + * @brief Class that defines the lbm module which is the interface between the 1D solver and OLB. +*/ +template +class lbmOocSimulator : public lbmMixingSimulator { + +using DESCRIPTOR = olb::descriptors::D2Q9<>; +using NoDynamics = olb::NoDynamics; +using BGKdynamics = olb::BGKdynamics; +using BounceBack = olb::BounceBack; + +using ADDESCRIPTOR = olb::descriptors::D2Q5; +using ADDynamics = olb::AdvectionDiffusionBGKdynamics; +using NoADDynamics = olb::NoDynamics; + +private: + + std::shared_ptr> tissue; + std::string organStlFile; ///< The STL file of the CFD domain. + + std::shared_ptr> organStlReader; + std::shared_ptr> organStl2Dindicator; + + T Vmax; + + void readOrganStl(); + + void prepareNsLattice(const T omega) override; + + void prepareAdLattice(const T omega, int speciesId); + +public: + /** + * @brief Constructor of an lbm module. + * @param[in] id Id of the module. + * @param[in] name Name of the module. + * @param[in] pos Absolute position of the module in _m_, from the bottom left corner of the microfluidic device. + * @param[in] size Size of the module in _m_. + * @param[in] nodes Map of nodes that are on the boundary of the module. + * @param[in] openings Map of the in-/outlets of the module. + * @param[in] stlFile STL file that describes the geometry of the CFD domain. + * @param[in] charPhysLength Characteristic physical length of the geometry of the module in _m_. + * @param[in] charPhysVelocity Characteristic physical velocity of the flow in the module in _m/s_. + * @param[in] alpha Relaxation factor for the iterative updates between the 1D and CFD solvers. + * @param[in] resolution Resolution of the CFD mesh in gridpoints per charPhysLength. + * @param[in] epsilon Convergence criterion for the pressure values at nodes on the boundary of the module. + * @param[in] relaxationTime Relaxation time tau for the LBM solver. + */ + lbmOocSimulator(int id, std::string name, std::string stlFile, std::shared_ptr> tissue, std::string organStlFile, std::shared_ptr> cfdModule, + std::unordered_map*> species, std::unordered_map> openings, ResistanceModel* resistanceModel, T charPhysLenth, T charPhysVelocity, + T alpha, T resolution, T epsilon, T relaxationTime=0.932, T adRelaxationTime=0.932); + + /** + * @brief Initialize an instance of the LBM solver for this module. + * @param[in] dynViscosity Dynamic viscosity of the simulated fluid in _kg / m s_. + * @param[in] density Density of the simulated fluid in _kg / m^3_. + */ + void lbmInit(T dynViscosity, T density) override; + + /** + * @brief Prepare the LBM geometry of this instance. + */ + void prepareGeometry() override; + + /** + * @brief Prepare the LBM lattice on the LBM geometry. + */ + void prepareLattice() override; + + /** + * @brief Write the vtk file with results of the CFD simulation to file system. + * @param[in] iT Iteration step. + */ + void writeVTK(int iT); + +}; + +} // namespace arch diff --git a/src/simulation/simulators/olbOoc.hh b/src/simulation/simulators/olbOoc.hh new file mode 100644 index 0000000..29a1262 --- /dev/null +++ b/src/simulation/simulators/olbOoc.hh @@ -0,0 +1,239 @@ +#include "olbOoc.h" +#include + +namespace sim{ + +template +lbmOocSimulator::lbmOocSimulator ( + int id_, std::string name_, std::string stlFile_, std::shared_ptr> tissue_, std::string organStlFile_, std::shared_ptr> cfdModule_, + std::unordered_map*> species_, std::unordered_map> openings_, ResistanceModel* resistanceModel_, T charPhysLength_, + T charPhysVelocity_, T alpha_, T resolution_, T epsilon_, T relaxationTime_, T adRelaxationTime_) : + lbmMixingSimulator(id_, name_, stlFile_, cfdModule_, species_, openings_, resistanceModel_, charPhysLength_, charPhysVelocity_, + alpha_, resolution_, epsilon_, relaxationTime_, adRelaxationTime_), + tissue(tissue_), organStlFile(organStlFile_) + { + this->cfdModule->setModuleTypeLbm(); + this->fluxWall.try_emplace(int(0), &this->zeroFlux); + } + +template +void lbmOocSimulator::lbmInit (T dynViscosity, T density) { + + this->setOutputDir(); + this->initValueContainers(); + this->initNsConverter(dynViscosity, density); + this->initAdConverters(density); + this->initNsConvergeTracker(); + this->initAdConvergenceTracker(); + + Vmax = (*tissue->getVmax(0))*this->getAdConverter(0).getPhysDeltaT(); + + #ifdef VERBOSE + std::cout << "[lbmSimulator] lbmInit " << this->name << "... OK" << std::endl; + #endif +} + +template +void lbmOocSimulator::prepareGeometry () { + + bool print = false; + + #ifdef VERBOSE + print = true; + #endif + + this->readGeometryStl(print); + this->readOpenings(); + readOrganStl(); + + this->geometry->clean(print); + this->geometry->checkForErrors(print); + + #ifdef VERBOSE + std::cout << "[lbmSimulator] prepare geometry " << this->name << "... OK" << std::endl; + #endif +} + +template +void lbmOocSimulator::prepareLattice () { + + /** + * Prepare the NS lattice + */ + const T omega = this->converter->getLatticeRelaxationFrequency(); + prepareNsLattice(omega); + + /** + * Prepare the AD lattices + */ + for (auto& [speciesId, converter] : this->adConverges) { + const T adOmega = this->getAdConverter(speciesId).getLatticeRelaxationFrequency(); + prepareAdLattice(adOmega, speciesId); + } + + /** + * Initialize the integral fluxes for the in- and outlets + */ + this->initPressureIntegralPlane(); + this->initFlowRateIntegralPlane(); + for (auto& [adKey, LatticeAD] : this->adLattices) { + this->initConcentrationIntegralPlane(adKey); + } + + /** + * Initialize the lattices + */ + this->initNsLattice(omega); + for (auto& [speciesId, converter] : this->adConverges) { + this->initAdLattice(speciesId); + } + + std::cout << "[lbmSimulator] prepare lattice " << this->name << "... OK" << std::endl; + + this->prepareCoupling(); + + std::cout << "[lbmSimulator] prepare coupling " << this->name << "... OK" << std::endl; +} + +template +void lbmOocSimulator::writeVTK (int iT) { + + bool print = false; + #ifdef VERBOSE + print = true; + #endif + + olb::SuperVTMwriter2D vtmWriter( this->name ); + // Writes geometry to file system + if (iT == 0) { + olb::SuperLatticeGeometry2D writeGeometry (this->getLattice(), this->getGeometry()); + vtmWriter.write(writeGeometry); + vtmWriter.createMasterFile(); + this->vtkFile = olb::singleton::directories().getVtkOutDir() + olb::createFileName( this->name ) + ".pvd"; + } + + if (iT % 1000 == 0) { + + olb::SuperLatticePhysVelocity2D velocity(this->getLattice(), this->getConverter()); + olb::SuperLatticePhysPressure2D pressure(this->getLattice(), this->getConverter()); + olb::SuperLatticeDensity2D latDensity(this->getLattice()); + vtmWriter.addFunctor(velocity); + vtmWriter.addFunctor(pressure); + vtmWriter.addFunctor(latDensity); + + vtmWriter.write(iT); + + // write all concentrations + for (auto& [speciesId, adLattice] : this->adLattices) { + olb::SuperLatticeDensity2D concentration( this->getAdLattice(speciesId) ); + concentration.getName() = "concentration " + std::to_string(speciesId); + vtmWriter.write(concentration, iT); + } + + // write vtk to file system + this->vtkFile = olb::singleton::directories().getVtkOutDir() + "data/" + olb::createFileName( this->name, iT ) + ".vtm"; + this->converge->takeValue(this->getLattice().getStatistics().getAverageEnergy(), print); + } + if (iT %1000 == 0) { + #ifdef VERBOSE + std::cout << "[writeVTK] " << this->name << " currently at timestep " << iT << std::endl; + #endif + } + + this->converge->takeValue(this->getLattice().getStatistics().getAverageEnergy(), print); + + if (iT%100 == 0) { + if (this->converge->hasConverged()) { + this->isConverged = true; + } + } + +} + +template +void lbmOocSimulator::readOrganStl () { + // Define Organ area + organStlReader = std::make_shared>(organStlFile, this->converter->getConversionFactorLength()); + organStl2Dindicator = std::make_shared>(*organStlReader); + this->geometry->rename(1, 3, *organStl2Dindicator); +} + +template +void lbmOocSimulator::prepareNsLattice (const T omega) { + + this->lattice = std::make_shared>(this->getGeometry()); + + // Initial conditions + std::vector velocity(T(0), T(0)); + olb::AnalyticalConst2D rhoF(1); + olb::AnalyticalConst2D uF(velocity); + + // Set lattice dynamics + this->lattice->template defineDynamics(this->getGeometry(), 0); + this->lattice->template defineDynamics(this->getGeometry(), 1); + this->lattice->template defineDynamics(this->getGeometry(), 2); + this->lattice->template defineDynamics(this->getGeometry(), 3); + + // Set initial conditions + this->lattice->defineRhoU(this->getGeometry(), 1, rhoF, uF); + this->lattice->iniEquilibrium(this->getGeometry(), 1, rhoF, uF); + + // Set lattice dynamics and initial condition for in- and outlets + for (auto& [key, Opening] : this->moduleOpenings) { + if (this->groundNodes.at(key)) { + setInterpolatedVelocityBoundary(this->getLattice(), omega, this->getGeometry(), key+3); + } else { + setInterpolatedPressureBoundary(this->getLattice(), omega, this->getGeometry(), key+3); + } + } +} + +template +void lbmOocSimulator::prepareAdLattice (const T adOmega, int speciesId) { + + std::shared_ptr> adLattice = std::make_shared>(this->getGeometry()); + + // Initial conditions + std::vector zeroVelocity(T(0), T(0)); + olb::AnalyticalConst2D rhoF(1); + olb::AnalyticalConst2D rhoZero(0); + olb::AnalyticalConst2D uZero(zeroVelocity); + + // Set AD lattice dynamics + adLattice->template defineDynamics(this->getGeometry(), 0); + adLattice->template defineDynamics(this->getGeometry(), 1); + adLattice->template defineDynamics(this->getGeometry(), 2); + adLattice->template defineDynamics(this->getGeometry(), 3); + + // Set initial conditions + adLattice->defineRhoU(this->getGeometry(), 0, rhoZero, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 0, rhoZero, uZero); + adLattice->defineRhoU(this->getGeometry(), 1, rhoF, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 1, rhoF, uZero); + adLattice->defineRhoU(this->getGeometry(), 2, rhoF, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 2, rhoF, uZero); + adLattice->defineRhoU(this->getGeometry(), 3, rhoF, uZero); + adLattice->iniEquilibrium(this->getGeometry(), 3, rhoF, uZero); + + adLattice->template defineField(this->getGeometry(), 0, uZero); + adLattice->template defineField(this->getGeometry(), 1, uZero); + adLattice->template defineField(this->getGeometry(), 2, uZero); + adLattice->template defineField(this->getGeometry(), 3, uZero); + + for (auto& [key, Opening] : this->moduleOpenings) { + adLattice->template defineDynamics(this->getGeometry(), key+3); + adLattice->iniEquilibrium(this->getGeometry(), key+3, rhoF, uZero); + adLattice->template defineField(this->getGeometry(), key+3, uZero); + } + + // Add wall boundary + olb::setFunctionalRegularizedHeatFluxBoundary(*adLattice, adOmega, this->getGeometry(), 2, this->fluxWall.at(0), this->fluxWall.at(0)); + + // Add Tissue boundary + olb::setFunctionalRegularizedHeatFluxBoundary(*adLattice, adOmega, this->getGeometry(), 3, &Vmax, tissue->getkM(speciesId)); + + this->adLattices.try_emplace(speciesId, adLattice); + +} + +} // namespace arch \ No newline at end of file