[jan]

Dear Rookie,

I am glad to hear that it works now and thank you very much for providing the working source code.

Best regards,
Jan

[jan]

Dear Rookie,

Loading the saved data will indeed restart the simulation at the first time step, but with the saved fluid data. If you want to continue from the stopped time step, I’d suggest to write this information additionally to another file and then read the time step from this file when loading the fluid data. Then use it to set iT accordingly. See [1] for an introduction to file input/output, if necessary.

Unfortunately, there is currently no implemented method for saving and loading the particle data. You could do it manually, for example by writing data (global particle ID, position, velocity, force, …) to files and loading them later, creating new particles and updating their data with the input from the file.

Best regards,
Jan

[1]: https://cplusplus.com/doc/tutorial/files/

[jan]

Dear avrachan,

If you only want to define the bottom wall, using IndicInverse<T, DESCRIPTOR::d> is wrong. This indicator inverts the indicator that is passed to it, i.e. thebottomboxyou define is the hole in a wall that takes up the rest of the indicator’s bounding box (defined byminandmax). Therefore, I suggest that you only use theIndicatorCuboid directly:

Vector<T,2> extend_mL(lengthX,L/2.);
Vector<T,2> origin_mL(0.,-L/2);
solidBoundaries.push_back( SolidBoundary<T, DESCRIPTOR::d>(
std::make_unique<IndicatorCuboid2D<T, DESCRIPTOR::d>>(extend_mL,origin_mL), 2, wallContactMaterial));

The error you saw is caused by the particles being completely immersed in the walls. Therefore, a contact is detected and a force is calculated. Full overlap is unrealistic, and the contact model is not valid in this case either, making the force obviously incorrect and likely causing unexpected behavior.

Regarding the second question. You’re absolutely right, you could wrap that part with another if condition that checks if the corresponding location is near a wall (e.g. using the circumference radius of the geometry). However, in most cases the signed distance functions should be fairly cheap to evaluate, and in my tests the contact detection on the lattice wasn’t the limiting factor, the force calculation was usually more expensive, but that might admittedly change if many complex walls are involved.

Best regards,
Jan

• This reply was modified 1 week, 1 day ago by jan.

[jan]

Dear Rookie,

1.
The mentioned code should save and then immediately load the saved fluid data. However, that seems unnecessary, because you could just keep using the data that is already stored in the SuperLattice. Instead, you could load the solutions if they exist and then skip the fluid calculations:

if ( !( superLattice.load(  “bstep3d.checkpoint”  ) ) ) {
  // Calculate fluid solution

  // And save afterwards
  sLattice.save( “bstep3d.checkpoint” );
}
else {
  // Do something with the saved fluid data
}

The if condition is important here because it ensures that the fluid calculations are only performed if the checkpoint data doesn’t exist.
In your case, if you’d like to always load the fluid solution and then do the same calculations again from that point, for example to divide long simulations into shorter segments, this should work:

sLattice.load(“bstep3d.checkpoint”);
// Calculate fluid solution (main loop)
// And save afterwards
sLattice.save( “bstep3d.checkpoint” );
// Terminate program and restart

2.
Of course, you can also do the fluid calculation at the same time. To do this you just have to call superLattice.collideAndStream(); again (and maybe also setBoundaryValues( superLattice, converter, iT, fluidMaxPhysT, superGeometry ); if necessary). Note, however, that without back-coupling the fluid velocity doesn’t really change anymore, so it seems unnecessary to keep calculating it. Only if you add some dynamic boundary conditions it might be necessary.

Best regards,
Jan

• This reply was modified 1 week, 2 days ago by jan.

[jan]

Dear jakschee,

as Adrian already mentioned, the function is correct.

I think the parameters are self-explanatory, but the indPtr would be a std::shared_ptr to a STLreader in your case. The latticeSpacing should be equal to or less than the Δx used in the simulation. It refers to the resolution of the second lattice used to derive particle information.

If you’re interested in the underlying method, please refer to the corresponding publications of our group [1-4].

The restriction to convex particles only applies to the discrete contact model [5]. However, for small particle volume fractions, an explicit contact model may not be necessary [6].

Best regards,
Jan

[1]: 10.1016/j.compfluid.2018.02.027
[2]: 10.3390/computation9020011
[3]: 10.3390/computation9040040
[4]: 10.5445/IR/1000136875
[5]: 10.1016/j.partic.2022.12.005
[6]: 10.1016/j.jocs.2024.102263

[jan]

Dear Jijo,

I believe it would be possible to add explicit functionality for this, but the corresponding indicator, smooth indicator and creator function are not currently implemented.

You could use `setResolvedArbitraryShape2D’ as you mentioned, but you should be aware that it discretizes the geometry and therefore introduces some errors, and it’s also likely to be a bit more computationally expensive.

The missing rotation matrix comes from the PARTICLETYPE. In the dkt2d example, typedef ResolvedCircleWithContact2D PARTICLETYPE; is set. As you can see, it explicitly mentions circle because a circle is rotation invariant. So if you want a different shape, please use typedef ResolvedParticleWithContact2D PARTICLETYPE; instead. See src/particles/descriptor/particleDescriptorAlias.h for more aliases. You could also just combine different fields you need in the same way as it’s done there.

Best regards,
Jan

[jan]

Dear Javad,

as Jijo pointed out, you could try to use HLBM by using the PorousParticleKupershtokhForcedBGKdynamics. Here you could use the force term as in the ForcedBGKdynamics in the multicomponent examples, while simultaneously including the particles.

For the particulate flow the closest starting point would be the SettlingCubes3d example. However, it should also work to start with one of the multicomponent examples.

Best regards,
Jan

[jan]

Dear Jijo,

if you solely need a force field for the coupling, then perhaps try PorousParticleKupershtokhForcedBGKdynamics and see if that works.

Best regards,
Jan

[jan]

Dear Jijo,

gravity is applied as a (constant) external acceleration, see accExt and apply_gravity in examples/particles/dkt2d. So the changes of the fluid density don’t affect the particle.

Also, currently the external acceleration is defined for the whole particle system. Updating it within the ParticleManager would therefore update it for all particles. In your particular case, I think it would be best to manually apply the gravity at each timestep. For example, you could evaluate the density at the current particle position if an averaged density is sufficient. If you need different forces on the surfaces, then a more sophisticated approach would be necessary.

Best regards,
Jan

[jan]

Dear Rookie,

if the Spring School doesn’t work for you, then maybe becoming a member of the consortium (https://www.openlb.net/consortium/) is an alternative option, because unfortunately I can’t provide such detailed support here.

Best regards,
Jan

[jan]

Dear Rookie,

Henn has used this code extensively, so I suggest you check out his dissertation (DOI: 10.5445/IR/1000061601).

For more in-depth discussions, you should consider attending the next occurrence of our annual [Spring School](https://www.openlb.net/spring-school-2024/).

Best regards,
Jan

[jan]

Dear Rookie,

I may have misunderstood you, but I suggest that you still read the files “one by one”. However, you can use a simple loop in c++ to do this for you by getting a list of files with the c++ function mentioned above.

But I want to emphasize that this is just a guess that it will work in your case, because I don’t have enough knowledge about your problem/case.

For more in-depth discussions, you should consider attending the next occurrence of our annual [Spring School](https://www.openlb.net/spring-school-2024/).

Best regards,
Jan

[jan]

Dear Rookie,

unfortunately, I’m not aware of what exactly you want to achieve and which functions you are currently using, but I’m glad to hear that you found a way to make it work. Though I don’t understand the issue with loading them separately as you could just use c++ features to iterate over all existing vtu files as mentioned here: https://stackoverflow.com/questions/612097/how-can-i-get-the-list-of-files-in-a-directory-using-c-or-c

Best regards,
Jan

[jan]

Dear Rookie,

I meant this part:

_minPhys(sg.getStatistics().getMinPhysR(2)),
_maxPhys(sg.getStatistics().getMaxPhysR(2)),

I haven’t checked the actual implementation, but usually this is used to evaluate where the fluid starts and ends, so that a particle leaving the fluid can be placed on the other side. However, usually 2 is for a wall, so this wouldn’t return the start and end of the fluid, but of the wall. This could (but doesn’t have to) cause your particles to behave wrong or “strangely” when they cross the periodic boundary.

Best regards,
Jan

[jan]

Dear Rookie,

removing the const keyword shouldn’t have any impact as long as you don’t try to manipulate the values that are stored in it.

Could you please clarify what you replaced with int overlap = 2; and why?

The overlap defines the size of the area in which the communication of data happens. For communication, the blocks must so that storage for the necessary neighborhood data is provided and accessible. Usually an overlap of 2 is sufficient.

Best regards,
Jan

[Author]

Posts