DEM in combination with other geometries

Applications
1

In the “WaveTank” test case, the beach is defined mathematically as a linear slope. Currently (Dec. 2019), I’m working with the “next” branch, which defines the geometry of the beach as a rotated “addBox” feature with a fixed boundary fill type. The test case runs as expected, but I’m now trying to modify the beach geometry.

The beach I’m trying to include is uniform along the y-axis but has an elevation defined by a high order polynomial function of the x-axis. In my initial approach, I was going to try approximating the beach function using a piecewise linear function with very small step sizes and build the beach as in “WaveTank”. However, this proved to be rather tedious and I quickly searched for an alternative approach.

The simplest alternative I could come up with was to build a DEM of the beach and include it following the “DEMExample” test case. Specifically, I tried to build the geometry following the “WaveTank” example then add a beach DEM at the end of the flume. However, this does not appear to work and I have a few comments and questions…

I tried to add the DEM, then shift its origin further down the wave flume with the following

GeometryID dem = addDEM(dem_file);
shift(dem,1.75,0,0);
GeometryID fluid_box = addDEMFluidBox(H,dem);

The “shift” causes program failure. Commenting it out initializes beach DEM as expected, only it is not located at the intended position. It’s as if “shift” does not apply to DEM geometries. I do see that “TopoCube” has a “shift” function, so I don’t know what is causing the failure.

More importantly, the inclusion of the DEM seems to unexpectedly override, or deactivate, other geometries. Again, I’m using the “next” branch, in which the “WaveTank” domain is initialized with a fluid box in the following way

GeometryID fluid = addBox(GT_FLUID, FT_SOLID, m_origin, 1.75, ly, H);

When a DEM is included in the problem, this fluid box is not added. I tried to make sure the boxes don’t overlap and it still fails to generate.

I then tried to make the entire domain a DEM, using “addDEMFluidBox”, however I could not add additional geometries to the domain, such as a cylinder. That being said, the wave maker is generated without any problems.

This leads me to believe that “addBox” and the DEM do work together, but the compatibility depends on the geometry type?

In summary, these are NOT compatible with the DEM approach

  • “addBox(GT_FLUID,…”
  • “addCylinder(GT_FIXED_BOUNDARY,…”

whereas this is compatible

  • “addBox(GT_MOVING_BODY,…”

Overall, I think this DEM stuff is WAY over complicated. Couldn’t the creation of a DEM be treated exactly as an “addBox”, where there is no conflicts with including other geometries? In simplest form, a DEM is nothing more than a collection of boxes of varying heights. Right now it seems to operate as an isolated method of defining the geometry with a code structure that has an added layer of complication.

I’ve tried numerous approaches, all of which have failed. Can someone offer advice on how to build a nonlinear beach in the “WaveTank” example?

2

After some long drawn out trial and error cases, here’s an update…

  1. The DEM does work with other features, such as “addBox(GT_FLUID,…” and “addCylinder(GT_FIXED_BOUNDARY,…”. However, it appears that they must be added AFTER calling
    GeometryID dem = addDEM(dem_file);
    shift(dem,1.75,0,0);
    GeometryID fluid_box = addDEMFluidBox(H,dem);
    It is as if there is an order dependency on how the geometry is constructed?

  2. The DEM “shift” does work (not sure what was going wrong prior), however there does appear to be a bug, or at least it doesn’t behave the way you would think it should. I’ve found that the DEM will repeat itself over the domain. In other words, if the domain is 2 meters long, and the DEM is 0.5 meters long, there will be 4 adjacent copies of the DEM over the domain. The “shift” appears to act as a “window” or a " boxcar function" which starts at the shift location and then spans the width of the DEM. Only what lies in the “window” will be generated. Continuing my example, if I define the shift at 0.25, then GPUSPH will fill the domain with the last half of the DEM file (0.25-0.5), then see the next “copy” and fill the remainder of the domain with the first half of the DEM file (0-0.25). In other words, the window is staggered over copies 1 and 2, which is easily confirmed by changing the shift value. I should note that there does appear to be some boundary around each copy which causes some loss of information.

  3. The observation above (2) does not however solve my problem, as my DEM width is not an integer factor of the domain length. Therefore, I will always get an incomplete version of the DEM data. As a solution, I added the DEM without the “shift”, then offset the domain origin at -1.75,0,0 (following the example from original post). This works perfect in my case, because I want the DEM at the end of the wave flume. I believe this works because “TopoCube” is hard coded to have the geometry origin at 0,0,0 (which I think think this should change and the origin should be an input like “addBox”). Basically, I left the DEM at 0,0,0 which also prevents any weird boundary errors, then added my objects and wavemaker in the negative portion of the domain.

Overall, I’ve been able to get the geometries in the correct position and they do appear to initialize correctly. However, I’m still facing one problem. The particles are not respecting the DEM as a boundary. When the simulation starts, they just fall to the floor as if they were initialized in mid air. I tried using “DEMExample” with the same DEM file and the particles at least respect the DEM. The first few time steps are weird though…the particles shoot all over and generate a lot of “noise” when it should just sit there as still water. It’s as if the particles are falling after initialization then bouncing off the boundaries. I think this was addressed in the new “WaveTank” example in the “next” branch, so maybe it’s something that also needs to be addressed in “addDEMFluidBox”? Not sure. Advice on how to fix that would be helpful.

Going back to the main problem though, I don’t know why the particles would fall out of mid air in one example vs the other. I started building my problem from the “WaveTank” example, which I stripped everything from to match the domain in the “DEMExample”. Problem still occurs, which leads me to believe it’s not a domain issue. I think it has to be a setting or something? The “WaveTank” and “DEMExample” are two different frameworks/numerics/physics so I’m not really sure where to look. The “DEMExample” has minimal setup, which might be favorable for debugging.

I’m really hoping some one can give some feedback on this last issue. I’m kind of stumped as to why those particles would fall out of mid air?

3

Starting to feel like I’m talking to myself here…

Anyway, an update for those who travel this path. It turns out in the “next” branch, a “DUMMY_BOUNDARY” was introduced as opposed to the “LJ_BOUNDARY”. I don’t know what the “DUMMY_BOUNDARY” does, why it was introduced, what it fixed, nothing. No comments, docs, nothing. That being said, the WaveTank.cu test case does look much improved from the “master” branch. What I found is that the DEM geometry does not like this new boundary. I submitted a bug fix and hopefully it gets resolved in the future. (search for “WaveMaker.cu” and “DEM” posts on some of the latest issues https://github.com/GPUSPH/gpusph/issues)

As far as what I’m doing now…I ended up having to flip the domain around and do everything backwards. For some reason, shifting the origin backwards was causing some issues. I don’t remember what off the top of my head, but I ended up not doing that. Instead, I kept the DEM at the origin, then moved the wave maker down the flume. This required that I change the wave maker normal and angle stuff, which was not too bad. I can give more feedback on that if need be. As I said before, I then changed the boundary type and everything seems to be working now.

4

Hello @GWAVE, sorry for the late reply and thanks for looking into this a lot on your own.

You seem to have hit several issues at once.

One is that the order of definition of geometries, that is important due to the current mechanism for intersections/particle erasure, that operate on the previously defined geometries only. FWIW, I’m not a big fan of this approach either, but until we get to the redesigned geometry library with first-class union/subtraction/intersection of geometries, this is what we have.

Secondly, concerning the non-LJ boundary models: the DUMMY boundary is a particle model (like DYN boundaries), but with improved equations for imposing boundary conditions. It shouldn’t affect the DEM logic per se.

The DEM can work in two different ways. One ways is “geometrically”: it doesn’t get filled with particles, but you need the ENABLE_DEM flag for the framework. This is only supported with the LJ_BOUNDARY model. The other way is by actually define boundary particles covering the DEM, which is done by giving it a FT_BORDER filltype. This should work regardless of the boundary model used (modulo bugs).

You also seem to have hit an issue with the shift for the DEM. This is most likely a bug that we need to address. I’ve noticed that you’ve opened new tickets on GitHub for this and other issues, thank you very much, this’ll help us keep track of it.

5

@giuseppe.bilotta Thanks for your help! One last thing on this subject…

I was following the logic presented in the “DEMExample.cu” where it appears to follow the “geometric” approach. The framework is set to “ENABLE_DEM” and the DEM is added using “addDEM”, “addDEMFluidBox”, etc. In this case, you are correct in saying that it is only supported with the “LJ_BOUNDARY” model. This was a “bug” that I reported because I personally don’t think the DEM concept should not be limited to a specific boundary.

That being said, I was originally working with the “DUMMY_BOUNDARY” (no DEM) and I did notice that the boundary is a bit smoother. I would like to continue to use this boundary model, but I also require a DEM. You said above that there is another way to implement the DEM? In my experience, I get critical errors if I don’t “ENABLE_DEM” while calling other DEM routines. Can you elaborate a bit more on how I can create the DEM using method 2 you mentioned? I know what “FT_BORDER” does, but I don’t know what routine to pass this to. I should reiterate that my DEM can not be built from primitives, hence why I’m trying to use this functionality. If I can get the DEM to work with the “DUMMY_BOUNDARY”, it would be great.

6

It’s theoretically possible to support other boundary models with DEM and planes. For dummy and dynamic, this would be using the ghost particle approach, with virtual particles mirroring the real particles across the planes and with properties set based on the expected boundary condition. As usual, the lack of the feature is mostly due to lack of manpower.

To switch to the particle representation for the DEM, you can specify the fill type as third argument for addDEM (the first argument is the file name, the second argument is optional and specifies the file type, the third argument is optional too and specifies the fill type). So you should be able to do something like:

addDEM("your-file-name", DEM_FMT_ASCII /* or _VTK or _XYZ */, FT_BORDER);
7

@giuseppe.bilotta I’m a bit confused on the current feature status…

(Just to reiterate, I’m using “next” branch)

Are you saying that it is currently possible to build DEM using the particle representation, or this is something that can be done, but is not implemented at the moment? (i.e. “it is on the TODO list”)

In the second paragraph of your response, you provide some instruction. This made me believe there is current support. Here are the parts of the code that I focused on …

SETUP_FRAMEWORK(
               viscosity<SPSVISC>,
                boundary<DUMMY_BOUNDARY>,
                add_flags<ENABLE_DEM | ENABLE_PLANES>
        ).select_options(
                RHODIFF,
                use_planes, add_flags<ENABLE_PLANES>()
        );

GeometryID dem = addDEM("topo_flip.txt",DEM_FMT_ASCII,FT_BORDER);
  1. I tried it with and without “ENABLE_DEM”
  2. I tried it with and without a call to “addDEMFluidBox”

I consistently keep getting an “error” which says:

Generating problem particles…
FATAL: TopoCube::FillIn not implemented !

I traced this back through the source code…

grep -nr "TopoCube::FillIn" ./src/*
./src/geometries/TopoCube.h:157:	{ throw std::runtime_error("TopoCube::FillIn not implemented !"); }

The source code around line 157 in the file above reads…

void FillIn(PointVect& points, const double dx, const int layers) override
{ throw std::runtime_error("TopoCube::FillIn not implemented !"); }

so it’s not really an “error” but more of a known lack of support. This leads me to believe that your second paragraph is not a set of instructions, but rather a pseudo-code on how it could be implemented?

Am I right in my understanding…this is possible but currently not implemented? Or am I missing something?

8

Ouch, I had forgotten that this was missing, which means that DUMMY or DYN aren’t supported at all with DEMs currently. On the upside, this should be easier to implement.

(BTW, when using the border filling, you should not add the ENABLE_DEM flag; that’s only needed when using the geometric approach.)

9

@giuseppe.bilotta

In the past, you’ve pointed to the fact that many things can be done in theory, but there is a lack of support to implement all features upon request. I totally get it and I commend the GPUSPH team on what they’ve provided to this point.

If you claim

then this might be an opportunity for me to contribute. Starting with something simple would be a good way for me to learn the code structure and from what I learn, it may enable me to contribute more in the future. If you could provide

  1. a general overview of what objectives need to be accomplished to meet the goal
  2. which files need to be addressed (so I’m not trying to dissect the entire code on my own)
  3. a possible feature within the code that mirrors the same feature (maybe the “box” geometry or something) so I can follow along
  4. any other relevant information for development (the more you share, the less I would have to ask)

If you provide this information, then I just might be able to carry out the task of implementing this feature. I will review what you provide, cross reference it with the source code, and assess if it is something I feel I can honestly accomplish. If so, I will take it on the challenge. If I think it is too difficult, to the point where I have to ask you questions every step of the way, then I won’t take up your time.

Point being, if it is as easy as you make it sound, it might be a good learning opportunity for you to share with me. We could also start a new topic here on the Discourse to document the process of implementing a new feature. I don’t know how many people are contributing to the code development, but I think this may be a task that is not trivial to the community?

Let me know what you think

10

Hello @GWAVE, I appreciate your commitment.

The geometry associated to the DEM is called a TopoCube (cube with topography), and it’s defined in the TopoCube.cc and TopoCube.h under src/geometries (all geometry files are there).

The FillIn method is currently defined in the .h file with a simple throw for the exception you’re seeing. For a proper implementation, you should remove the current stub from the .h file and put a proper definition in the .cc file. You can base the implementation of FillIn on the current definition of the FillDem method (around line 475 of the .cc file).

Currently, what FillDem does is to define an (x, y) grid of points with spacing dx and then for each (x, y) coordinate find the DEM height and define a point there. For FillIn, you want to use the same logic, but instead of a single point at each (x, y) point you want to add layers point, starting from the DEM-computed z and moving down by dx for each. A simple for loop with each generated Point at z - i*dx should work.

Let me know if you need any more support.