diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9813ad0b4b023adf5c3100f712f2a8044430c1dc..507228da27f30ff986bebfe855fa6ce0a16b5984 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 ab5d625edd2baf3c101b6f19d944bd98e1458095..5b4d0c4def46dbe9fe795d106735aabb3cd4c9f2 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 ab9c48124ddcf87d9248744a0006e6d6477b42ad..a475ecaf03d0dec0dfe49075bc45aa01f1cc22c1 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 1850e63ec0d79ac565aed1c2a4598f4876f18a0e..db0009894b69897ee00ec7f4ef22cbdc145bbbd5 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 30d4f07082efd08f158bcb36e16322af7f7231e2..b5eb33610f41a82f60d76f2dea3e50c90fb17771 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>;