Adding Complex Boundaries

[ryan-bournes]

This update adds complex boundaries to the diffusion grid. This will allow the user to specify the boundary type separately on all 6 sides of the simulation space (Neumann, Dirichlet or Periodic).

This can be done by changing the diffusion_boundary_condition in the simulation parameters to 'Complex' and then setting the new simulation parameters diffusion_bc_x_min, diffusion_bc_x_max, diffusion_bc_y_min, diffusion_bc_y_max, diffusion_bc_z_min and diffusion_bc_z_max to either Neumann, Dirichlet or Periodic.

This update also comes with two new functions to the diffusion grid class: SetBoundaryCondition_neumann and SetBoundaryCondition_dirichlet. If the user is using complex boundaries and are using both Neumann and Dirichlet boundaries in the same simulation, they can set the boundary condition for these types separately. For example: one can set x-min boundary to Dirichlet and x_max boundary to Neumann, then use these functions to set the Dirichlet boundary to 10 and Neumann boundary to 0 (closed).

Some important notes about this update:

• The Dirichlet boundary works slightly differently in complex boundaries compared to the normal Dirichlet boundaries. An illustration of this can be found below: instead of the voxel at the boundary being set to the boundary condition, the voxel just outside the boundary is instead. This still functions the same way normal Dirichlet boundaries do, albeit it means the simulation space will fill up with concentration slightly slower, since the boundary condition is located slightly further away.

  Normal Dirichlet: BC=S c_-1 | S c_1 c_2 .....| Complex Dirichlet: BC=S S | c_0 c_1 c_2 .....|

• The SetBoundaryCondition_neumann and SetBoundaryCondition_dirichlet functions are only used for Complex. The base SetBoundaryCondition is still used to define the boundary condition for the base Dirichlet and Neumann boundaries.

• When setting the boundary conditions in complex boundaries, there is fatal check if the value inputted for any of the boundaries is something other than Neumann. Dirichlet or Periodic.

As this feature provides quite a bit of flexibility to the user, it does come with some unit tests. I have designed four unit tests for the complex boundaries:

1. A death test to make sure initialising the grid fails when using complex boundaries and one of the values for the boundaries is unrecognised.

2. This test is functionally the same as the EulerDirichletBoundaries test except the concentration only comes from the top instead of all sides. x_min is set to Dirichlet = 1 and x_max is set to Neumann = 0 (closed). All other boundaries are set to Periodic. Check to make sure for the first few iterations, the concentration at the Neumann boundary is less than the Dirichlet boundary. Then after many iterations, check to see if the average concentration in the whole simulation space is the same as the Dirichlet boundary.

3. This test is functionally the same as the EulerPeriodicBoundary test, except we should only see the concentration of the source at y_min and y_max observe as Periodic boundaries. x_min and x_max are set to Neumann=0, y_min and y_max are set to Periodic, z_min and z_max are set to Dirichlet=0. A source is placed at each boundary, and the concentration next to the source inside simulation space and at the opposite side of the simulation are compared. The source at x_min, x_max, z_min, z_max should see a greater concentration within the space compared to opposite, while the source at y_min and y_max should have equal concentrations either side.

4. This test sets up a basic steady state environment, with concentration coming in at the Neumann boundary and going out at the Dirichlet boundary. We should see the diffusion gradient set up between these two boundaries overtime. y_min is set to Neumann = -1 and y_max is set to Dirichlet = 0. All other boundaries are set to Periodic. Check to make sure for the first few iterations, the total concentration in the simulation is increasing. Then after many iterations,check to see if the change in concentration is very small (<0.0001). For the whole simulation, check to see if the concentration next to the Neumann boundary is greater than the concentration at the Dirichlet boundary.

@TobiasDuswald and @nicogno I am requesting for both of you to be reviewers since you both are most familiar with the diffusion grid. There may be some areas of this code that can be trimmed down or simplified. If you have any questions about any part of this update then please let me know.

[TobiasDuswald]

Seems as if the unit tests fail on macOS, and the code is not formatted yet @ryan-bournes .

Generally, if we introduce these complex boundaries, I think this should be the default case, and all others should build up on it. That would conveniently lead to all unit tests testing the complex boundaries for different cases. For example, the Neumann case should just set 6 Neumann boundaries for complex boundaries, and so forth, rather than having their own stencil implementation.

Also, looking at the code, it's starting to get messy with plenty of if statements for different cases. Can we cast the update into a matrix form y_{n+1} = A y_n? We'd build (assemble) the sparse matrix A once for the provided coefficients and boundaries and then simply execute a sparse matrix-vector product for each step. This would also allow us to use intelMKL, for instance, to accelerate and optimize that part of the code further in the future. Any thoughts on that?

[nicogno]

Hi @ryan-bournes , I agree with @TobiasDuswald regarding the use of the complex boundaries as the default case.

As for the matrix form for the state update, that's a great idea as well! We can discuss further about that implementation.

[ryan-bournes]

Hi Both,

Thank you for the feedback. I've always struggled with writing code in the most optimal way so getting some ideas to improve the code layout is very helpful. I agree that designing the code so that it reads from a matrix of boundary values would be be much better than a bunch of IF statements.

As for making complex boundaries to be the default case, I agree that it would be much better but I am not sure if that could affect any backwards compatibility? As each boundary type is its own separate function, the diffusion grid class has DiffuseWithDirichlet and DiffuseWithNeumann as well as Diffuse. By default, the continuum operation, all the demos and all the unit tests use Diffuse but I am not sure if any other projects use DiffuseWithDirichlet or the other functions directly?

Currently when running Diffuse, it simply locates which function to use based on what boundary type the parameter is set to, and this could be changed so the running Diffuse always uses the complex boundary function. When the user sets the boundary type to Neumann, it just sets all 6 sides to Neumann and the complex function should behave the same as the Neumann function. I still think it would be worth keeping the complex unit tests as they specifically test if the boundary can still set different types to the 6 sides of the environment.

In the meantime, I can investigate the issue this code has on mac. I use linux so I am not sure what the exact problem is but it may have something to do with the clamp function.

[TobiasDuswald]

Matrix

I may have phrased that misleadingly; I wondered if we could encode the stencil update together with the boundaries into a large, sparse matrix A to update the discretized concentration values u with a sparse matrix-vector multiplication.

Interface / compatibility

I can not exclude it, but theoretically, nobody should call DiffuseWithXYZ. I think we do it somewhere in the unit tests for convenience. Maybe these DiffuseWithXYZ should even become private in the long run. But I meant that all other boundaries should fall back to the "complex" boundaries, i.e., not change the interface, but have them call the complex boundaries while treating all boundaries as the method suggests.

clamp

I suppose you use std::clamp, this is part of the C++ standard and therefore highly unlikely to explain differences between linux / macos.

[sonarcloud]

Quality Gate failed

Failed conditions
5.8% Duplication on New Code (required ≤ 3%)

See analysis details on SonarCloud