From c12c66a581d68c918e97266b2995d7ee8ec4a6aa Mon Sep 17 00:00:00 2001
From: Lukas Riedel <mail@lukasriedel.com>
Date: Tue, 30 Jun 2020 23:18:15 +0000
Subject: [PATCH] Use BCGS_AMG_SSOR linear solver for FV computations

This selects the BCGS_AMG_SSOR linear solver backend for FV computations
in both the Richards and the Transport model.
---
 CHANGELOG.md                            |  5 +++-
 dune/dorie/model/richards/richards.cc   | 17 ++++++++++---
 dune/dorie/model/richards/richards.hh   |  6 +++--
 dune/dorie/model/transport/transport.cc | 33 +++++++++++++++++++------
 dune/dorie/model/transport/transport.hh | 28 +++++++++++++++------
 5 files changed, 67 insertions(+), 22 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9813ad0b..507228da 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,10 @@
 
 ## Unreleased
 
-<!-- Add new entries here! -->
+### Changed
+
+* Linerar solver for finite volumes changed from `AMG_4_DG` to
+    `BCGS_AMG_SSOR` !204
 
 ## 2.0.0 (2020-05-14)
 
diff --git a/dune/dorie/model/richards/richards.cc b/dune/dorie/model/richards/richards.cc
index ab5d625e..5b4d0c4d 100644
--- a/dune/dorie/model/richards/richards.cc
+++ b/dune/dorie/model/richards/richards.cc
@@ -154,9 +154,20 @@ void ModelRichards<Traits>::operator_setup()
 	igo = std::make_unique<IGO>(*go0,*go1);
 
 	// --- Solvers ---
-	lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
-	lscc = std::make_unique<LSCC>();
-	ls = std::make_unique<LS>(*igo,*cc,*lsgfs,*lscc,1000,0,true,true);
+
+	// FV linear solver
+	if constexpr (order == 0) {
+		// Solver args: GFS, maxiter, verbose, reuse, superlu
+		ls = std::make_unique<LS>(*gfs, 1000, 0, true, true);
+	}
+	// DG linear solver
+	else {
+		lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
+		lscc = std::make_unique<LSCC>();
+		// Solver args: GO, constraints, LSGFS, LSGFS constraints, maxiter, verbose,
+		//              reuse, superlu
+		ls = std::make_unique<LS>(*igo, *cc, *lsgfs, *lscc, 1000, 0, true, true);
+	}
 
 	pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
 	pdesolver->setParameters(inifile.sub("NewtonParameters"));
diff --git a/dune/dorie/model/richards/richards.hh b/dune/dorie/model/richards/richards.hh
index ab9c4812..a475ecaf 100644
--- a/dune/dorie/model/richards/richards.hh
+++ b/dune/dorie/model/richards/richards.hh
@@ -130,8 +130,10 @@ struct ModelRichardsTraits : public BaseTraits
 	/// Solution vector type
 	using U = typename IGO::Traits::Domain;
 	/// Linear solver types
-	using LS = Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<IGO,CC,LSGFS,LSCC,
-		Dune::PDELab::CG2DGProlongation,Dune::SeqSSOR,Dune::BiCGSTABSolver>;
+	using LS =
+		std::conditional_t<order == 0,
+			Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<IGO>,
+			Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<IGO, CC, LSGFS, LSCC, Dune::PDELab::CG2DGProlongation, Dune::SeqSSOR, Dune::BiCGSTABSolver>>;
 	/// Non-linear solver types
 	using PDESOLVER = Dune::PDELab::Newton<IGO,LS,U>;
 	/// Time stepping scheme
diff --git a/dune/dorie/model/transport/transport.cc b/dune/dorie/model/transport/transport.cc
index 1850e63e..db000989 100644
--- a/dune/dorie/model/transport/transport.cc
+++ b/dune/dorie/model/transport/transport.cc
@@ -182,16 +182,33 @@ void ModelTransport<Traits>::operator_setup()
     explicit_igo = std::make_unique<ExplicitIGO>(*go0,*go1);
 
   // --- Solvers ---
-  lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
-  lscc = std::make_unique<LSCC>();
+  // Initialize helper spaces for DG linear solver only
+  if constexpr (order > 0)
+  {
+    lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
+    lscc = std::make_unique<LSCC>();
+  }
+
   if (ts_param->implicit())
-    implicit_ls = std::make_unique<ImplicitLS>(
-      *implicit_igo, *cc, *lsgfs, *lscc, 1000, 0, true, true
-    );
+  {
+    if constexpr (order == 0)
+      implicit_ls = std::make_unique<ImplicitLS>(*gfs, 1000, 0, true, true);
+    else
+    {
+      implicit_ls = std::make_unique<ImplicitLS>(
+        *implicit_igo, *cc, *lsgfs, *lscc, 1000, 0, true, true);
+    }
+  }
   else
-    explicit_ls = std::make_unique<ExplicitLS>(
-      *explicit_igo, *cc, *lsgfs, *lscc, 1000, 0, true, true
-    );
+  {
+    if constexpr (order == 0)
+      explicit_ls = std::make_unique<ExplicitLS>(*gfs, 1000, 0, true, true);
+    else
+    {
+      explicit_ls = std::make_unique<ExplicitLS>(
+        *explicit_igo, *cc, *lsgfs, *lscc, 1000, 0, true, true);
+    }
+  }
 
   // --- Time Step Operators ---
   if (ts_param->implicit()){
diff --git a/dune/dorie/model/transport/transport.hh b/dune/dorie/model/transport/transport.hh
index 30d4f070..b5eb3361 100644
--- a/dune/dorie/model/transport/transport.hh
+++ b/dune/dorie/model/transport/transport.hh
@@ -158,14 +158,26 @@ struct ModelTransportTraits : public BaseTraits
   /// Solution vector type
   using U = typename ImplicitIGO::Traits::Domain;
   /// Linear solver type
-  using ImplicitLS = Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<
-    ImplicitIGO, CC, LSGFS, LSCC, Dune::PDELab::CG2DGProlongation,
-    Dune::SeqSSOR, Dune::BiCGSTABSolver
-  >;
-  using ExplicitLS = Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<
-    ExplicitIGO, CC, LSGFS, LSCC, Dune::PDELab::CG2DGProlongation,
-    Dune::SeqSSOR, Dune::BiCGSTABSolver
-  >;
+  using ImplicitLS = std::conditional_t<
+    order == 0,
+    Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<ImplicitIGO>,
+    Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<ImplicitIGO,
+                                            CC,
+                                            LSGFS,
+                                            LSCC,
+                                            Dune::PDELab::CG2DGProlongation,
+                                            Dune::SeqSSOR,
+                                            Dune::BiCGSTABSolver>>;
+  using ExplicitLS = std::conditional_t<
+    order == 0,
+    Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<ExplicitIGO>,
+    Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<ExplicitIGO,
+                                            CC,
+                                            LSGFS,
+                                            LSCC,
+                                            Dune::PDELab::CG2DGProlongation,
+                                            Dune::SeqSSOR,
+                                            Dune::BiCGSTABSolver>>;
   /// Methods computing the time step
   using ImplicitSLPS
     = Dune::PDELab::StationaryLinearProblemSolver<ImplicitIGO, ImplicitLS, U>;
-- 
GitLab