Add first unit test

[mrp089]

Problem

When implementing new features, having unit tests would be very useful. This is especially true for solid material models since they have a clearly defined interface:

• in: deformation gradient $\mathbf F$
• out: 2PK stress $\mathbf S$, tangent matrix $\mathbb{C}$

It's also straight-forward to come up with analytical solutions for materials to test against.

Solution

We will add a single unit test, starting with the Neo-Hookean material under different loading conditions. Once we have a framework integrated into GitHub actions, we can expand it to other currently implemented hyperelastic materials (since we need to cover them anyway). Finally, we can use that framework to add new material models (#2, #33, #49).

Additional context

@yuecheng-yu will add a .cpp and a cmake file. @mrp089 will integrate the executable into GitHub actions. @MatteoSalvador, input welcome!

Code of Conduct

• [X] I agree to follow this project's Code of Conduct and Contributing Guidelines

[ktbolt]

@mrp089 As I recall I said let's not do this because it is a lot of work and the code was not implemented with unit testing in mind.

There are lots of other more important things that need to be done.

If @yuecheng-yu has some time then he can work on https://github.com/SimVascular/svFSIplus/issues/30 or update our material models with those of svFSI.

[mrp089]

Reopening this after talking with @ktbolt, so we have a place to discuss and exchange ideas.

My PhD lab is using GoogleTest (or gtest) in their (closed-source) multi-physics FEM code. Here is their 30 min video guide how to include GoogleTest (slides). This was part of a Summer of Testing series. Here is a quick start on including it using CMake.

A challenge in the current svFSI+ will be creating mock objects, especially for ComMod.

[mrp089]

@yuecheng-yu has prototyped a first unit test using GoogleTest and calling get_pk2cc on a solid material model. He's currently cleaning it up after running into issues with cmake and relative path definitions. The first application could be testing for a stress-free reference configuration in a couple of material models.

[yuecheng-yu]

@mrp089 I think I find the "secrete" location of Google Test. Finally, I don't choose to separate the unit test's CMakeLists, because it never works at any directories (as we discussed today). Instead, the main program's CMakeLists (.../svFSI/CMakeLists.txt) extends to build both Google Test and Unit Test. During compilation, after the test cases (hello_test.cpp, test.cpp) built, Google Test will "pre-run" all of them from the directory: build/svFSI-build/Source/svFSI. This will be the current_path for all relative path. After make, instead of run the executable directly (this is what I made mistake I think), run ctest from build/svFSI-build/Source/svFSI and Google Test will automatically run all test cases. All the test results or print output will be generated at build/svFSI-build/Source/svFSI/Testing/Temporary/LastTest.log. In this way, I can ensure the consistency of current_path at both compilation time and run time. The test source codes and input files are separated at tests/unitTests.

[yuecheng-yu]

I tried to replace gtest_discover_tests(test_exe) by add_test(NAME TestFamily COMMAND test_exe) in svFSI's CMakeLists (.../svFSI/CMakeLists.txt). It can avoid "pre-run" during make, but remain the function to use ctest to run test cases.

[mrp089]

Google Mock might be a way to get around reading an input file to create a (Simulation and) ComMod object. We could mock everything in ComMod that's not actually needed by get_pk2cc. Is this what you had in mind, @ktbolt?

[ktbolt]

@mrp089 The idea is to make the unit tests as simple and independent of other objects as possible, don't want to need to modify a bunch of unit tests when a class changes. And you don't want to read in files.

Creating mock objects is a way to implement this, most unit test frameworks support creating mock objects.

[aabrown100-git]

I've started expanding the unit testing for material models a bit. You can find the latest commits here: https://github.com/aabrown100-git/svFSIplus/tree/material_unit_tests

I wrote tests to check the consistency of the Mooney-Rivlin 2nd PK stress $\mathbf{S}$ and material elasticity tensor $\mathbb{C}$ produced by get_pk2cc() for a random deformation gradient $\mathbf{F}$ using finite differences. First, I wrote a function to compute the strain energy density function $\Psi$ as a function of $\mathbf{F}$ (Psi_MR()). Then, I'm checking that the following equalities are true

TEST 1: see TEST(TestMR, S_3D)

$\mathbf{S}(\mathbf{F}):\Delta\mathbf{E} == \Delta\Psi$

where $\Delta\mathbf{E} = \mathbf{E}(\mathbf{F} + \Delta\mathbf{F}) - \mathbf{E}(\mathbf{F})$ and $\Delta\Psi =\Psi(\mathbf{F} + \Delta\mathbf{F}) -\Psi(\mathbf{F})$ and $\Delta\mathbf{F}$ is a random perturbation of the deformation gradient.

(To be clear, $\mathbf{S}(\mathbf{F})$ is obtained from get_pk2cc() and $\Psi(\mathbf{F})$ is provided by the user).

Checking this for many random $\Delta\mathbf{F}$ should be enough to verify that $\mathbf{S}(\mathbf{F})$ is correct.

TEST 2: see TEST(TestMR, CC_3D)

$\mathbb{C}(\mathbf{F}):\Delta\mathbf{E} == \Delta\mathbf{S}$

where $\Delta\mathbf{E} = \mathbf{E}(\mathbf{F} + \Delta\mathbf{F}) - \mathbf{E}(\mathbf{F})$ and $\Delta\mathbf{S} =\mathbf{S}(\mathbf{F} + \Delta\mathbf{F}) -\mathbf{S}(\mathbf{F})$ and $\Delta\mathbf{F}$ is a random perturbation of the deformation gradient.

(To be clear, $\mathbb{C}(\mathbf{F})$ is obtained from get_pk2cc() and $\mathbf{S}(\mathbf{F})$ is also obtained from get_pk2cc(). Thus, we're only checking whether, within get_pk2cc(), $\mathbb{C}(\mathbf{F})$ is consistent with $\mathbf{S}(\mathbf{F})$ )

Checking this for many random $\Delta\mathbf{F}$ should be enough to verify that $\mathbb{C}(\mathbf{F})$ is correct.

TEST 3: see TEST(TestMR, S_3D_direct)

$\mathbf{S}(\mathbf{F}) == \hat{\mathbf{S}}(\mathbf{F})$

where $\mathbf{S}(\mathbf{F})$ is obtained from get_pk2cc() and $\hat{\mathbf{S}}(\mathbf{F})$ is obtained by finite differencing. Specifically $\hat{\mathbf{S}}(\mathbf{F}) = \mathbf{F}^{-1}\hat{\mathbf{P}}(\mathbf{F})$ where $\hat{P}_{iJ} = \frac{\Psi(\mathbf{F}^{(iJ)}) - \Psi(\mathbf{F})}{\Delta}$ and $\mathbf{F}^{(iJ)}$ is $\mathbf{F}$ with its $iJ$ component perturbed by $\Delta$.

(I tried to do a similar thing to compute $\mathbb{C}(\mathbf{F})$ using finite differencing, but ran into some issues. This test is probably redundant with TEST 1, but I'm keeping it for now since it seems to be working.)

I think all of these tests are useful because they only require the user to implement $\Psi(\mathbf{F})$, which is usually explicitly stated in papers, whereas $\mathbf{S}(\mathbf{F})$ and $\mathbb{C}(\mathbf{F})$ typically must be derived by hand, which can lead to errors.

I could use some help thinking through the best way to organize these tests. For example, we currently have a single class called UnitTestIso which is meant to be a class to test isotropic material models. I think it would be nice to have a generic class for testing any hyperelastic material model (including volumetric strain energy functions), and for each specific hyperelastic material have a derived class which contains the specific strain energy function for that material. It should also contain the correct parameters and matType so that it obtains the appropriate $\mathbf{S}$ and $\mathbb{C}$ from get_pk2cc(). Then, I think we should be able to write these 3 tests for any material model in a very concise manner.

[mrp089]

I think that's an excellent idea, @aabrown100-git! I'd be happy to help make this a pull request.

I have two ideas:

1. Can you move the FD function from the tests to the solid material models? This way, we could define an input flag to use the FD tangent for material evaluation. This would enable quickly prototyping materials. On the test side, we would compare the two function outputs, analytic and FD.
2. In my experience, it can sometimes be hard to say at what tolerance an FD tangent is right or wrong (for a given epsilon). A more reliable but slightly more complicated way is to compare the error for at least two epsilons. You could then check the rate of convergence of your FD scheme (e.g., linear with forward/backward FD, quadratic with central FD) instead of the absolute error.

Let me know what you think!

[ktbolt]

@aabrown100-git A good start.

Some comments:

• Since @lpapamanolis is working on a new class-based material model would it be better to wait on this? It will need to be totally redone with a new implementation.

• Test names should clearly describe what the test is testing. TEST(TestMR, S_3D_direct) should be named something like TEST(TestMooneyRivlenMateral, StressTensor3D_Valid), not really sure what S is.TestMooneyRivlenMateral might should be TestMaterialModels_MooneyRivlin.

• Use a test fixture to create data used by several tests.

• Unit tests should not be very complex, a few lines of code so it is clear what is being tested. In the case of testing get_pk2cc<3>() you would just set F, nFn, fN,ya_g

  
  Array F{{1.0, 0.0, 0.0},
        {0.0, 1.0, 0.0},
        {0.0, 0.0, 1.0}};
  
  nFn = 1;
  
  Array Fn{{1.0, 1.0, 0.0},
         {1.0, 1.0, 0.0},
         {0.0, 1.0, 1.0}};
  
  double ya_g = 1.0;
  
  mat_models_carray::get_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, fN, ya_g, S, Dm);

and then compare S and Dm to what you would expect their values to be given the values for F, nFn, fN,ya_g.

You can set F, nFn, fN,ya_g to test certain edge cases and test error handling.

- The unit test files should be organized like the source code they are testing.

[aabrown100-git]

Thanks @mrp089 and @ktbolt. I've made a few changes based on your suggestions. Currently, the code is a bit disorganized, but I think it will become cleaner after we make some decisions.

@mrp089 Regarding your two ideas

1. Currently, I'm only calculating the PK2 stress uses finite differences in the strain energy function. I'm running into some issues doing a similar thing for the material elasticity tensor. Maybe we can meet sometime to discuss?
2. Yes, this is a good idea. I think I'll come back to this when I've organized the rest of the code.

@ktbolt

1. Yes, it would be better to wait for @lpapamanolis, but I'm hoping to implement and verify a new material model soon.
2. I made test names more descriptive following your suggestion
3. I added test fixtures for each material model
4. I consolidated some code to make units tests more concise, but there's more work to be done

[aabrown100-git]

I've reorganized the material model tests and added a variety of tests for the Holzapfel-Ogden material model. On my computer, these tests pass, but they sometimes fail due to randomness in generating test deformation gradient tensors. Following @mrp089 suggestion, I will think about how to test the rate of convergence of the finite difference approximations instead.

I implemented these tests primarily to verify the Holzapfel-Ogden material model implementation for a cardiac mechanics benchmark problem I want to solve with svFSIplus. I think the tests are currently sufficient for that.

@ktbolt @mrp089 I think merging it with svFSIplus would be useful, but we may also want to wait until @lpapamanolis finishes the material model restructuring. If we want to merge it soon, please look over the code when you get a chance and let me know what can be improved.

[ktbolt]

@aabrown100-git You should be using the C++ random library and not rand(). The C++ random library has better algorithms and produces the same results on all platforms, rand() could produce different results on different platforms.

If the test are useful then I say go ahead and merge unless this breaks something.

[aabrown100-git]

@ktbolt Thanks for the advice Dave, I just modified the code to switch to the C++ random library.

I'll submit a merge and we can continue the conversation there!

[aabrown100-git]

I've added some more tests (testing ustruct material model implementations and testing volumetric penalty model implementations for both struct and ustruct). I'll submit a new merge request.

[aabrown100-git]

I'm expanding the tests to test the order of convergence of the finite difference estimates for S and CC. You can find the latest here. One issue I'm dealing with is that as delta decreases, the error decreases as expected, but eventually it bottoms out around 1e-6~1e-7 (due to round-off error?). I am computing the order of convergence by finding the slope of log(error) vs. log(delta), so if this bottomed-out error is included the order of convergence is computed incorrectly. Is there some way to ignore this bottomed-out error when calculating order of convergence? Also, any ideas why the error bottoms out so high (in double precision shouldn't it bottom out around 1e-16?).

[ktbolt]

@aabrown100-git It's not clear to me what the motivation is for adding more unit tests. With 50 open svFSIplus Issues I have to think that there are more important things to work on than adding another unit test that will in future need to be refactored.

But if you think this is important then create a new Issue for it so it can be more clearly tracked.

[aabrown100-git]

@ktbolt I'm interested in improving the unit tests to have confidence in the material model implementations, especially because I and others want to run some cardiac mechanics benchmark cases in svFSIplus.

I'll create a new issue for this.

[aabrown100-git]

I think perhaps we can close this issue?

[ktbolt]

@aabrown100-git Gaining confidence in the material models is definitely a good thing, have at it!