For my purposes, I need to be able to set a certain inflow profile, for example a variant of a Poiseuille flow where the velocity changes at every time step.
Currently, this is possible when using generated sweeps, but only with a few nasty hacks, for both making it temporally varying and spatially varying. To be able to use a DynamicUBB, which was generated via
ubb_dynamic = UBB(lambda *args: None, dim=stencil.D)
ubb_data_handler = UBBAdditionalDataHandler(stencil, ubb_dynamic)
# UBB with user-defined velocity profile
generate_boundary(ctx, "DynamicUBB", ubb_dynamic, lbm_method,
additional_data_handler=ubb_data_handler, target=target)
one needs to pass a std::function< Vector3< real_t >(const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&)
to the constructor of the DynamicUBB type. The easiest way to make one is via a functor, for which a template looks like this:
class InflowProfile
{
public:
Vector3< real_t > operator()( const Cell& pos, const shared_ptr< StructuredBlockForest >& SbF, IBlock& block ) {
// return velocity vector depending on the cell location in the SbF
return getVelocityVector(pos, SbF, block);
}
};
Temporally varying inflow profile
The problem is that, in the current version of waLBerla/lbmpy, the additional data handler is only being called once before the running of the simulation when the generated member function
template<typename FlagField_T>
void DynamicUBB::fillFromFlagField( const shared_ptr<StructuredBlockForest> & blocks, ConstBlockDataID flagFieldID,
FlagUID boundaryFlagUID, FlagUID domainFlagUID)
is being called. After these values have been set, there is currently no easy way to update them at every time step. As a workaround, one can add the aforementioned member function as a addFuncAfterTimeStep
to the SweepTimeLoop and add a TimeTracker
object to the functor, which allows to update the values depending on the current time step. This approach does appear to incur a performance penalty, as the member function DynamicUBB::fillFromFlagField
appears to do some work which should be unnecessary after an initial run of it. I think, but haven't tested it, that this also makes running on the GPU not possible or significantly slowed down.
Spatially varying inflow profile
Currently, the cells being passed to the functors operator()
are different from the ones I expected to get. The problem is in the following generated code:
// for every cell in the block
if ( isFlagSet( it.neighbor(1, 0, 0 , 0 ), boundaryFlag ) )
{
auto element = IndexInfo(it.x(), it.y(), it.z(), 0 );
Vector3<real_t> InitialisatonAdditionalData = elementInitaliser(Cell(it.x(), it.y(), it.z()), blocks, *block);
element.vel_0 = InitialisatonAdditionalData[0];
element.vel_1 = InitialisatonAdditionalData[1];
element.vel_2 = InitialisatonAdditionalData[2];
// snip
}
// Similar code for all directions, e.g. 27 in total for a D3Q27 stencil
Where the function elementInitialiser
is the our functors operator()
.
In essence, what it being checked if the current cell has a neighbor, in the case of a inflow boundary, for which the inflowFlag is set. If that is the case, apply the functor to this cell, meaning the cell of which a neighboring cell has the inflowFlag set, not the cell with the inflowFlag self, and store the results (i.e. the computed velocity vector) in the appropriate location.
Currently, I am working around this problem by manually adding the direction to the cell on which the functor is applied. Unfortunately, due to the generated nature of the code, this is not a real solution...
In my use case, certain noise values are generated for every cell for every time step in the inflow boundary via a Python program which are then being read from a .json file. These values are then put into a std::unordered_map<GlobalCellCoordinatesTuple, Vector3>
and then read during the boundary treatment. Without the aforementioned fix there is a miss match between the cell coordinates I put in the std::unordered_map<..,..>
(the ones for which the inflowFlag is being set) and the cell coordinates that get passed to the functor, resulting in the wrong behavior.
Suggestions
TimeTracker
to the functor and then can update the values in an easy manner.operator()
to be 1) the neighboring cells as is done right now, or 2) they are the actual cells for which the flag which is checked is set.Dear Markus,
Thanks for the extensive reply, it definitely makes sense to do it this way.
I think for me the Blender approach would be better at this stage, as it would allow me to focus on quickly prototyping some implementations. I fear that writing a custom vertexFaceToColor
function for every mesh will become quickly intractable when dealing with more or complexer meshes.
Or is there an easy way to paint the appropriate vertices? I had a look through the OpenMesh documentation about traversing the mesh and couldn't figure out how to do so in a nice way.
If the author is still available and has time, then that would be great! From what I read in the thesis is that his changes are twofold: A Blender plug-in to also export colors to the .obj file and a way to generate a certain configuration file. If I were to have access to these additions, I think I should be able to come a long way from there!
And no worries, always happy to give back a bit of my time! Small amount of effort for a big gain for the community.
(I assume that this is the correct place to ask question about using waLBerla. If not, apologies!)
I'm a master student in Mathematics at TUM, writing my master thesis on using LBM to simulate the blood flow near cerebral aneurysms, and intend on using waLBerla as my LBM library, as I really like what I've seen so far.
First question beforehand: What is nowadays the preferred way to use waLBerla? I've read in some other gitlab issues that using lbmpy and pystencils are usually the way to go, either to generate C++ code or run 'natively' in python, but didn't see any documentation on the website on when/how to make this choice.
To do so, I have some .obj file with a triangle mesh of an artery near a cerebral aneurysm. When I simply use this object instead of, for example, the bunny and make some other small changes (i.e. exclude the scaling of the domain as done in the tutorial) and compile/run the program, I get the warning message that my mesh is not watertight and that is has boundary edges. As my mesh indeed is not watertight, I adjusted that in Blender, see also the attached pictures. The warning remains, so it has to do with the boundary edges.
My simulation should be relatively straightforward, I have one inflow, two outflows and on the geometry a noslip boundary condition. I understand that I must let the library know where these boundaries are, but I can't find how to do this in the documentation and all of the example applications don't do similar things.
I have found that the 2014 master thesis by Al Hasnine and the 2013 waLBerla paper about blood flow simulations on the coronary artery tree are very similar to my topics, but was unable to find the code corresponding to these simulations.
Are there any examples available on how to do such things? Or what would be the best strategy?
As I assume more people wanting to experiment with waLBerla may have similar questions, so I have no problem working the answer to this question out into a (small) tutorial!
Many thanks in advance!
Hi!
I would like to experiment a bit with waLBerla, but can't for the life of me figure out how to properly link against it. I asked for help on SO already, where someone mentioned that waLBerla does not provide a real CMake package config file.
More details can be found in the SO post about my setup so far, but as I'm pretty sure I'm doing things suboptimally, I was wondering what you could recommend. In general a tutorial would be nice and would definitely lower the barrier for new(er) users to experiment with waLBerla. If I manage to get it to work, I have no problem writing it up and make a merge request. Or include it in the setup instructions.
In short, say I want to compile the following code (from the first tutorial):
#include "core/Environment.h"
#include "blockforest/Initialization.h"
using namespace walberla;
int main( int argc, char ** argv )
{
Environment env( argc, argv );
// create blocks
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid(
uint_c( 3), uint_c(2), uint_c( 4), // number of blocks in x,y,z direction
uint_c(10), uint_c(8), uint_c(12), // how many cells per block (x,y,z)
real_c(0.5), // dx: length of one cell in physical coordinates
false, // one block per process? - "false" means all blocks to one process
false, false, false ); // no periodicity
return 0;
}
What would be the best possible project layout?
Thanks!
Since it works, I will close this issue. I have also added an answer to my SO post, so in the future people can find it a bit more easily!
Yes, I can see that file and it works, great!
When running tutorial 4, I got linker errors, because I didn't specify the mesh
module as a dependency. Do I always have to manually figure out which dependencies I have or is there a smarter way to do that?
Hi!
I would like to experiment a bit with waLBerla, but can't for the life of me figure out how to properly link against it. I asked for help on SO already, where someone mentioned that waLBerla does not provide a real CMake package config file.
More details can be found in the SO post about my setup so far, but as I'm pretty sure I'm doing things suboptimally, I was wondering what you could recommend. In general a tutorial would be nice and would definitely lower the barrier for new(er) users to experiment with waLBerla. If I manage to get it to work, I have no problem writing it up and make a merge request. Or include it in the setup instructions.
In short, say I want to compile the following code (from the first tutorial):
#include "core/Environment.h"
#include "blockforest/Initialization.h"
using namespace walberla;
int main( int argc, char ** argv )
{
Environment env( argc, argv );
// create blocks
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid(
uint_c( 3), uint_c(2), uint_c( 4), // number of blocks in x,y,z direction
uint_c(10), uint_c(8), uint_c(12), // how many cells per block (x,y,z)
real_c(0.5), // dx: length of one cell in physical coordinates
false, // one block per process? - "false" means all blocks to one process
false, false, false ); // no periodicity
return 0;
}
What would be the best possible project layout?
Thanks!