Exercise Model (DuMuX course)

The aim of this exercise is it to learn how to set up a new model (new system of equations). As an example, we implement a nonlinear diffusion equation model and apply it for image denoising. In this exercise, we will only consider the bare minimum of classes to successfully assemble and solve such a problem with DuMux. We also implement the model for a specific discretization method (the Box method: a vertex-centered finite volume method also known as control-volume finite element method with piece-wise linear basis functions).

Equation and Problem set-up

We consider a nonlinear diffusion equation on a domain \Omega \subset \mathbb{R}^2 and time interval (0,T] with boundary \partial\Omega given by

\begin{aligned}
\frac{\partial c}{\partial t} - \nabla \cdot (D(c) \nabla c) &= 0    && \mathrm{in}\; \Omega \times (0,T] \\
                         -D(c) \nabla c \cdot \boldsymbol{n} &= 0  && \mathrm{on}\; \partial\Omega \times (0,T] \\
                                                        c &= c^0  && \mathrm{for}\; t = 0,
\end{aligned}

where for a position \mathbf{x} \in \Omega and time t, c(\mathbf{x},t) denotes the unknown evolving concentration field and D(c) the solution-dependent diffusion coefficient.

In this exercise we want to model image denoising such that the primary variable c corresponds to the image intensity and D is given (due to Perona & Malik), as

D(c) = \frac{1}{1 + \left(\frac{|| \nabla c ||}{K}\right)^2},

with conductance K.

Task 1: Getting familiar with the diffusion example

When implementing a new model it is advisable to built upon an already existing model that might be similar and only needs modifications. Since the above equation is a nonlinear diffusion equation, we can built upon the diffusion model implemented in the diffusion example. The diffusion example also derives the discrete equations using the Box method as spatial discretization scheme.

Task 2: Set up the model

➡️ Copy the model.hh file from the diffusion example into dumux-course/exercises/exercise-model and choose appropriate class names.

Do also not forget to change the include guards at the beginning of header files (DUMUX_EXAMPLES_DIFFUSION_MODEL_HH). Include guards have to be unique in the entire application that you compile. (Otherwise some code will not be included.) An alternative is using #pragma once which is widely supported but not specified by the C++ standard.

First, the goal is to get the simulation running and then add improvements. For this, it is important to have a compiling test such that new changes can continuously be tested. Thus, in the computeFlux(...) function of the local residual, you can start with a hard-coded diffusion coefficient of 1.0 (linear diffusion coefficient function). (Replace problem.diffusionCoefficient() by 1.0 because our problem class in main.cc does not have a diffusionCoefficient() interface.)

Each model also needs to define a model type tag for setting model-specific properties.

➡️ Rename the one of the diffusion model (struct DiffusionModel {};) to NonlinearDiffusionModel. and use this new type tag for specializing properties, i.e. within the properties namespace

namespace Dumux::Properties {

// Defined properties, replace TTag::DiffusionModel by TTag::NonlinearDiffusionModel

} // end namespace Dumux::Properties

Task 3: Set up the test case

➡️ Open the file main.cc.

There you can find

1. The problem setup (implements the homogeneous Neumann boundary conditions (zero flux))
2. Test case specific property definitions
3. A function for setting the initial intensity values by using the noisy image data
4. The main program defining all steps of the program

For reading the image data (given as pgm file), the NetPBMReader reader is used

// Read the image
const auto imageFileName = getParam<std::string>("ImageFile");
const auto imageData = NetPBMReader::readPGM(imageFileName);

3.1: Make test type tag inherit properties from model type tag

➡️ Include the header model.hh
➡️ To use the new model for this test case, make the test type tag inherit properties from the model type tag

namespace Dumux::Properties::TTag {

struct NonlinearDiffusionTest
{
    // TODO: Set new model by replacing ModelTypeTag
    using InheritsFrom = std::tuple</*ModelTypeTag,*/BoxModel>;

    // Further property definitions
};
} // end namespace Dumux::Properties::TTag

3.2: Uncomment the commented code in main file

Without a model the previous code could not be compiled. The parts of the main function that need a working model are commented.

➡️ Remove the code comments to enable all simulation steps.

3.3: Run your application

We hope everything is implemented correctly.

➡️ Compile and run the test case by executing

cd build-cmake/exercises/exercise-model
make exercise_nonlineardiffusion
./exercise_nonlineardiffusion

➡️ Fix any compiler errors and warnings.

Task 4: Nonlinear diffusion

In the local residual class, instead of using a constant 1.0, we now want to implement the nonlinear diffusion function specified above. With a Dune::FieldVector you can get the 2-norm of the vector via v.two_norm() and its square with v.two_norm2(). (The reason for having two different interfaces is that the second one is more efficient.) Compute the diffusion coefficient and pass it to vmtv instead of 1.0. You can get the conductance K via the problem interface conductance.

➡️ Add the nonlinear diffusion coefficient computation to the computeFlux function

➡️ Compare the main.cc with the main.cc of the diffusion example. Notice that we use a Newton solver as we want to solve nonlinear equations.

As the DuMux assembler uses numeric differentiation to approximate the Jacobian matrix, you do not have to implement the derivatives of your residual by hand and implement them. This greatly simplifies implementing nonlinear equations.

Task 5: Running the image denosing test case

➡️ Compile and run like above.

With the conductance parameter Problem.Conductance (see params.input) you can regulate the preservation of edges. For large conductance, you get almost linear diffusion (Gaussian blurring). For small conductances some edges are preserved. The conductance acts as a threshold for the image intensity gradient detected as an edge feature.

The final result should look like this:

Task 6: Custom volume variables

The volume variables represent variables of the control volume. So far we have used the BasicVolumeVariables class that only stores the primary variables at control volumes. We have used the interface priVar to access the primary variables. In more complex models, we often want to

• assign names to our variables to make the local residual more readable
• evaluate constitutive laws and compute secondary variables from primary variables and use these in the local residual

In this task we will implement a custom volume variable class, while in this simple case it will only be for the first reason.

➡️ Implement your own ImageDenoisingVolumeVariables class (either within a new header file or directly in model.hh). As basis you can use the BasicVolumeVariables, see dumux/common/volumevariables.hh. Copy the class into model.hh and rename it.

Add a function imageIntensity() that returns the first primary variable.

➡️ In order to use the new ImageDenoisingVolumeVariables class in your model, set in your model-specific properties (see dumux-course/exercises/exercise-properties):

template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::NonlinearDiffusionModel>
{
    struct Traits
    {
        using PrimaryVariables
            = GetPropType<TypeTag, Properties::PrimaryVariables>;
    };
    // TODO: set type to your class
    using type = ...;
};

You can now use the interface imageIntensity in your local residual instead of the generic priVar interface.

In order to simplify the implementation of your custom volume variables, you can also inherit from BasicVolumeVariables and only implement the additional interface.

Another common modification is to add named indices to a model. To this end, we typically add an Indices struct to the ModelTraits in which the primary variable and equation indices are named (e.g. static constexpr int imageIntensityIdx = 0;). The model traits can be passed to the volume variables class via the Traits class.

➡️ Have a look at the solution of the exercise to see how this is usually implemented.

Names indices allow to address the entries of the local residual vector or the primary variables vector. You will see names indices in every DuMux model.