diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index ea2132c05712327a0923d05dee86e8bc9a32dd08..148f2c4ca79ae435782d3c4a690832a497a184cd 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/localgeodesicfeadolcstiffness.hh
@@ -103,7 +103,12 @@ energy(const typename Basis::LocalView& localView,
       localASolution[i] = aRaw[i];  // may contain a projection onto M -- needs to be done in adouble
     }
 
-    energy = localEnergy_->energy(localView,localASolution);
+    try {
+        energy = localEnergy_->energy(localView,localASolution);
+    } catch (Dune::Exception &e) {
+        trace_off(rank);
+        throw e;
+    }
 
     energy >>= pureEnergy;
 
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 0ccc13069583735c60c0f7af370ddc8dc77bca31..e70e6649f60d03a8e3e69ce99285c9812fbaad00 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -577,45 +577,54 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
 
             for (size_t j=0; j<newIterate.size(); j++)
                 newIterate[j] = TargetSpace::exp(newIterate[j], corr[j]);
-
-            energy  = assembler_->computeEnergy(newIterate);
-            energy = grid_->comm().sum(energy);
-
-            // compute the model decrease
-            // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
-            // Note that rhs = -g
-            CorrectionType tmp(corr.size());
-            tmp = 0;
-            hessianMatrix_->umv(corr, tmp);
-            modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-            modelDecrease = grid_->comm().sum(modelDecrease);
-
-            double relativeModelDecrease = modelDecrease / std::fabs(energy);
-
-            if (this->verbosity_ == NumProc::FULL and rank==0) {
-                std::cout << "Absolute model decrease: " << modelDecrease
-                          << ",  functional decrease: " << oldEnergy - energy << std::endl;
-                std::cout << "Relative model decrease: " << relativeModelDecrease
-                          << ",  functional decrease: " << (oldEnergy - energy)/energy << std::endl;
-            }
-            assert(modelDecrease >= 0);
-
-
-            if (energy >= oldEnergy and rank==0) {
-                if (this->verbosity_ == NumProc::FULL)
-                    printf("Richtung ist keine Abstiegsrichtung!\n");
+            try {
+                energy  = assembler_->computeEnergy(newIterate);
+            } catch (Dune::Exception &e) {
+                std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
+                std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
+                newIterate = x_;
+                solved = false;
+                energy = oldEnergy;
             }
-
-            if (energy >= oldEnergy &&
-                (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
-                if (this->verbosity_ == NumProc::FULL and rank==0)
-                    std::cout << "Suspecting rounding problems" << std::endl;
-
-                if (this->verbosity_ != NumProc::QUIET and rank==0)
-                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
-
-                x_ = newIterate;
-                break;
+            if (solved) {
+                energy = grid_->comm().sum(energy);
+
+                // compute the model decrease
+                // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+                // Note that rhs = -g
+                CorrectionType tmp(corr.size());
+                tmp = 0;
+                hessianMatrix_->umv(corr, tmp);
+                modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+                modelDecrease = grid_->comm().sum(modelDecrease);
+
+                double relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+                if (this->verbosity_ == NumProc::FULL and rank==0) {
+                    std::cout << "Absolute model decrease: " << modelDecrease
+                              << ",  functional decrease: " << oldEnergy - energy << std::endl;
+                    std::cout << "Relative model decrease: " << relativeModelDecrease
+                              << ",  functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+                }
+                assert(modelDecrease >= 0);
+
+
+                if (energy >= oldEnergy and rank==0) {
+                    if (this->verbosity_ == NumProc::FULL)
+                        printf("Richtung ist keine Abstiegsrichtung!\n");
+                }
+
+                if (energy >= oldEnergy &&
+                    (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
+                    if (this->verbosity_ == NumProc::FULL and rank==0)
+                        std::cout << "Suspecting rounding problems" << std::endl;
+
+                    if (this->verbosity_ != NumProc::QUIET and rank==0)
+                        std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+                    x_ = newIterate;
+                    break;
+                }
             }
         }
 
@@ -633,8 +642,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
 
             recomputeGradientHessian = true;
 
-        } else if (solved && (oldEnergy-energy) / modelDecrease > 0.01
-                    || std::abs(oldEnergy-energy) < 1e-12) {
+        } else if (solved && ((oldEnergy-energy) / modelDecrease > 0.01
+                    || std::abs(oldEnergy-energy) < 1e-12)) {
             // successful iteration
             x_ = newIterate;