diff -ruN exercises/exercise-model/CMakeLists.txt exercises/solution/exercise-model/CMakeLists.txt
--- exercises/exercise-model/CMakeLists.txt 2025-03-03 15:56:59.908625000 +0100
+++ exercises/solution/exercise-model/CMakeLists.txt 2025-03-03 15:56:59.911624948 +0100
@@ -3,5 +3,5 @@
dune_symlink_to_source_files(FILES images params.input)
-dumux_add_test(NAME exercise_nonlineardiffusion
+dumux_add_test(NAME exercise_nonlineardiffusion_sol
SOURCES main.cc)
diff -ruN exercises/exercise-model/main.cc exercises/solution/exercise-model/main.cc
--- exercises/exercise-model/main.cc 2025-03-03 15:56:59.908625000 +0100
+++ exercises/solution/exercise-model/main.cc 2025-03-03 15:56:59.911624948 +0100
@@ -29,6 +29,8 @@
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
+#include "model.hh"
+
/*!
* \ingroup NonlinearDiffusion
* \brief Test problem for image denoising using a nonlinear diffusion model
@@ -78,7 +80,7 @@
struct NonlinearDiffusionTest
{
// TODO: We need to set our model by replacing ModelTypeTag
- using InheritsFrom = std::tuple</*ModelTypeTag,*/BoxModel>;
+ using InheritsFrom = std::tuple<NonlinearDiffusionModel, BoxModel>;
using Scalar = double;
using Grid = Dune::YaspGrid<2>;
@@ -136,63 +138,62 @@
// the problem for the boundary conditions, a solution vector and a grid variables instance.
auto gridGeometry = std::make_shared<GridGeometry>(gridManager.grid().leafGridView());
- // TODO: Uncomment when the model is implemented
- // using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- // using Problem = GetPropType<TypeTag, Properties::Problem>;
- // using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
- // using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
-
- // auto problem = std::make_shared<Problem>(gridGeometry);
- // auto sol = createInitialSolution<SolutionVector>(*gridGeometry, imageData);
- // auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
- // gridVariables->init(sol);
-
- // // We initialize the VTK output module and write out the initial concentration field
- // VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
- // vtkWriter.addVolumeVariable([](const auto& vv){ return vv.priVar(0); }, "imageIntensity");
- // vtkWriter.write(0.0);
-
- // // We instantiate time loop using start and end time as well as
- // // the time step size from the parameter tree (`params.input`)
- // auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(
- // getParam<Scalar>("TimeLoop.TStart", 0.0),
- // getParam<Scalar>("TimeLoop.Dt"),
- // getParam<Scalar>("TimeLoop.TEnd")
- // );
-
- // // Next, we choose the type of assembler, linear solver and nonlinear solver
- // // and construct instances of these classes.
- // using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
- // using LinearSolver = SSORCGIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>;
- // using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
-
- // auto oldSol = sol; // copy the vector to store state of previous time step
- // auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol);
- // auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
- // Solver solver(assembler, linearSolver);
-
- // // time loop
- // timeLoop->start(); do
- // {
- // // assemble & solve
- // solver.solve(sol);
-
- // // make the new solution the old solution
- // oldSol = sol;
- // gridVariables->advanceTimeStep();
-
- // // advance to the time loop to the next step
- // timeLoop->advanceTimeStep();
-
- // // write VTK output (writes out the concentration field)
- // vtkWriter.write(timeLoop->time());
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+
+ auto problem = std::make_shared<Problem>(gridGeometry);
+ auto sol = createInitialSolution<SolutionVector>(*gridGeometry, imageData);
+ auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+ gridVariables->init(sol);
+
+ // We initialize the VTK output module and write out the initial concentration field
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
+ vtkWriter.addVolumeVariable([](const auto& vv){ return vv.imageIntensity(); }, "imageIntensity");
+ vtkWriter.write(0.0);
+
+ // We instantiate time loop using start and end time as well as
+ // the time step size from the parameter tree (`params.input`)
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(
+ getParam<Scalar>("TimeLoop.TStart", 0.0),
+ getParam<Scalar>("TimeLoop.Dt"),
+ getParam<Scalar>("TimeLoop.TEnd")
+ );
+
+ // Next, we choose the type of assembler, linear solver and nonlinear solver
+ // and construct instances of these classes.
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ using LinearSolver = SSORCGIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>;
+ using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+
+ auto oldSol = sol; // copy the vector to store state of previous time step
+ auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol);
+ auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
+ Solver solver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // assemble & solve
+ solver.solve(sol);
+
+ // make the new solution the old solution
+ oldSol = sol;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write VTK output (writes out the concentration field)
+ vtkWriter.write(timeLoop->time());
- // // report statistics of this time step
- // timeLoop->reportTimeStep();
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
- // } while (!timeLoop->finished());
+ } while (!timeLoop->finished());
- // timeLoop->finalize(gridGeometry->gridView().comm());
+ timeLoop->finalize(gridGeometry->gridView().comm());
return 0;
}// end main
diff -ruN exercises/exercise-model/model.hh exercises/solution/exercise-model/model.hh
--- exercises/exercise-model/model.hh 1970-01-01 01:00:00.000000000 +0100
+++ exercises/solution/exercise-model/model.hh 2025-03-03 15:56:59.911624948 +0100
@@ -0,0 +1,210 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux-Course contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+
+#ifndef DUMUX_EXERCISE_NONLINEAR_DIFFUSION_MODEL_HH
+#define DUMUX_EXERCISE_NONLINEAR_DIFFUSION_MODEL_HH
+
+// In the file `model.hh`, we define the volume variables, the model equations, and
+// set all default model properties.
+
+#include <dune/common/fvector.hh>
+#include <dumux/common/math.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/volumevariables.hh>
+#include <dumux/discretization/method.hh>
+
+// The property tag is simply an empty struct with the name `NonlinearDiffusionModel`
+namespace Dumux::Properties::TTag {
+//! The nonlinear diffusion model tag that we can specialize properties for
+struct NonlinearDiffusionModel {};
+} // end namespace Dumux::Properties::TTag
+
+namespace Dumux {
+
+/*!
+ * \ingroup NonlinearDiffusionModel
+ * \brief Volume averaged quantities required by the nonlinear diffusion model.
+ * This contains the quantities which are related to a control volume
+ */
+template <class Traits>
+class ImageDenoisingVolumeVariables
+{
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+public:
+ //! export the type used for the primary variables
+ using PrimaryVariables = typename Traits::PrimaryVariables;
+ using Indices = typename Traits::ModelTraits::Indices;
+
+ template<class ElementSolution, class Problem, class Element, class SubControlVolume>
+ void update(const ElementSolution& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ priVars_ = elemSol[scv.indexInElement()];
+ }
+
+ /*!
+ * \brief Returns the image intensity of the control volume.
+ */
+ Scalar imageIntensity() const
+ { return priVars_[Indices::imageIntensityIdx]; }
+
+ Scalar priVar(const int pvIdx) const
+ { return priVars_[pvIdx]; }
+
+ const PrimaryVariables& priVars() const
+ { return priVars_; }
+
+ Scalar extrusionFactor() const
+ { return 1.0; }
+
+private:
+ PrimaryVariables priVars_;
+};
+
+/*!
+ * \ingroup NonlinearDiffusionModel
+ * \brief Element-wise calculation of the residual and its derivatives
+ * for a nonlinear diffusion problem.
+ */
+template<class TypeTag>
+class NonlinearDiffusionModelLocalResidual
+: public GetPropType<TypeTag, Properties::BaseLocalResidual>
+{
+ using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables;
+ using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
+ using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using Indices = typename ModelTraits::Indices;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+public:
+ using ParentType::ParentType;
+
+ NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ NumEqVector storage;
+ storage[Indices::diffusionEqIdx] = volVars.imageIntensity();
+ return storage;
+ }
+
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
+ {
+ static_assert(DiscretizationMethods::isCVFE<typename GridGeometry::DiscretizationMethod>,
+ "This local residual is hard-coded to control-volume finite element schemes");
+
+ // Compute ∇c at the integration point of this sub control volume face.
+ const auto& fluxVarCache = elemFluxVarsCache[scvf];
+ Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ // v.axpy(a, w) means v <- v + a*w
+ gradConcentration.axpy(
+ volVars.imageIntensity(),
+ fluxVarCache.gradN(scv.indexInElement())
+ );
+ }
+
+ NumEqVector flux;
+
+ const auto K = problem.conductance();
+ const auto D = 1.0/(1.0 + gradConcentration.two_norm2()/(K*K));
+
+ // Compute the flux with `vtmv` (vector transposed times matrix times vector) or -n^T D ∇c A.
+ // The diffusion coefficient comes from the `problem` (see Part 2 of the example).
+ flux[Indices::diffusionEqIdx] = -1.0*vtmv(
+ scvf.unitOuterNormal(), D, gradConcentration
+ )*scvf.area();
+
+ return flux;
+ }
+};
+} // end namespace Dumux
+
+// The model properties
+namespace Dumux::Properties {
+
+// The type of the local residual is the class defined above
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::NonlinearDiffusionModel>
+{ using type = NonlinearDiffusionModelLocalResidual<TypeTag>; };
+
+// The default scalar type is double
+// we compute with double precision floating point numbers
+template<class TypeTag>
+struct Scalar<TypeTag, TTag::NonlinearDiffusionModel>
+{ using type = double; };
+
+// The model traits specify some information about our equation system.
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::NonlinearDiffusionModel>
+{
+ struct type
+ {
+ struct Indices
+ {
+ static constexpr int imageIntensityIdx = 0;
+ static constexpr int diffusionEqIdx = 0;
+ };
+
+ static constexpr int numEq() { return 1; }
+ };
+};
+
+// The primary variable vector has entries of type `Scalar` and is
+// as large as the number of equations
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::NonlinearDiffusionModel>
+{
+ using type = Dune::FieldVector<
+ GetPropType<TypeTag, Properties::Scalar>,
+ GetPropType<TypeTag, Properties::ModelTraits>::numEq()
+ >;
+};
+
+// The volume variables that are used in our model
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::NonlinearDiffusionModel>
+{
+ struct Traits
+ {
+ using ModelTraits
+ = GetPropType<TypeTag, Properties::ModelTraits>;
+ using PrimaryVariables
+ = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ };
+
+ using type = ImageDenoisingVolumeVariables<Traits>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff -ruN exercises/exercise-model/README.md exercises/solution/exercise-model/README.md
--- exercises/exercise-model/README.md 2025-03-03 15:23:33.245567666 +0100
+++ exercises/solution/exercise-model/README.md 1970-01-01 01:00:00.000000000 +0100
@@ -1,195 +0,0 @@
-# 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$, with boundary $\partial\Omega$, and time interval $(0,T]$ given by
-
-```math
-\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, $D(c)$ the solution-dependent diffusion coefficient, and $c^0$ the initial concentration field.
-
-In this exercise we want to model image denoising such that the primary variable $c$ corresponds to the **image intensity**, $c^0$ to the noisy image data, and $D$ is given (due to Perona & Malik), as
-```math
-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](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/tree/master/examples/diffusion).
-The diffusion example also derives the discrete equations using the Box method as spatial discretization scheme.
-
-## Task 2: Set up the model
-
-:arrow_right: 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](https://en.wikipedia.org/wiki/Include_guard)
-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`](https://en.wikipedia.org/wiki/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.
-
-:arrow_right: 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
-
-```c++
-namespace Dumux::Properties::TTag {
-
-// Defined properties, replace DiffusionModel by NonlinearDiffusionModel
-
-} // end namespace Dumux::Properties::TTag
-```
-
-## Task 3: Set up the test case
-
-:arrow_right: 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
-```c++
-// 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
-
-:arrow_right: Include the header `model.hh` <br>
-:arrow_right: To use the new model for this test case, make
-the `NonlinearDiffusionTest` type tag inherit properties from the model type tag
-
-```c++
-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.
-
-:arrow_right: Remove the code comments to enable all simulation steps.
-
-### 3.3: Run your application
-
-We hope everything is implemented correctly.
-
-:arrow_right: Compile and run the test case by executing
-
-```bash
-cd build-cmake/exercises/exercise-model
-make exercise_nonlineardiffusion
-./exercise_nonlineardiffusion
-```
-
-:arrow_right: Fix any compiler errors and warnings.
-
-## Task 4: Nonlinear diffusion
-
-In the local residual class, instead of using a constant diffusion coefficient of `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`.
-
-:arrow_right: Add the nonlinear diffusion coefficient computation to the `computeFlux` function
-
-:arrow_right: 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
-
-:arrow_right: 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:__
-<figure>
- <center>
- <img src="../../slides/img/exercise_model_mri_denoise.gif" alt="denoising"/>
- <figcaption> <b> Fig.1 </b> - Denosing of MRI image using nonlinear diffusion model.</figcaption>
- </center>
-</figure>
-
-## Task 6: Custom volume variables
-
-The volume variables represent variables related to control volumes. 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.
-
-:arrow_right: 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.
-
-:arrow_right: Add a function `imageIntensity()` that returns the first primary variable.
-
-:arrow_right: In order to use the new `ImageDenoisingVolumeVariables` class in your model, set in your model-specific properties (see `dumux-course/exercises/exercise-properties`):
-```c++
-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.
-
-:arrow_right: 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.