From 62e4d72cfc1073a841ee4c02b277c71bb075f1a5 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Wed, 4 Oct 2023 11:52:18 -0400
Subject: [PATCH] replaced nonlinear solver by linear solver

---
 festim/concentration/theta.py | 31 ++++++++++++++++++++-----------
 1 file changed, 20 insertions(+), 11 deletions(-)

diff --git a/festim/concentration/theta.py b/festim/concentration/theta.py
index aaa302122..11505eb72 100644
--- a/festim/concentration/theta.py
+++ b/festim/concentration/theta.py
@@ -88,26 +88,35 @@ def post_processing_solution_to_concentration(self):
         The attribute post_processing_solution is fenics.Product (if self.S is
         festim.ArheniusCoeff)
         """
-        du = f.TrialFunction(self.post_processing_solution.function_space())
-        J = f.derivative(self.form_post_processing, self.post_processing_solution, du)
-        problem = f.NonlinearVariationalProblem(
-            self.form_post_processing, self.post_processing_solution, [], J
+        problem = f.LinearVariationalProblem(
+            a=f.lhs(self.form_post_processing),
+            L=f.rhs(self.form_post_processing),
+            u=self.post_processing_solution,
+            bcs=[],
         )
-        solver = f.NonlinearVariationalSolver(problem)
-        # TODO these prms should be the same as in Simulation.settings I think
-        solver.parameters["newton_solver"]["absolute_tolerance"] = 1e-10
-        solver.parameters["newton_solver"]["relative_tolerance"] = 1e-10
-        solver.parameters["newton_solver"]["maximum_iterations"] = 50
+        solver = f.LinearVariationalSolver(problem)
         solver.solve()
 
     def create_form_post_processing(self, V, materials, dx):
+        """Creates a variational formulation for c = theta * S or theta**2 * S
+
+        Args:
+            V (fenics.FunctionSpace): the DG1 function space of the concentration field
+            materials (festim.Materials): the materials
+            dx (fenics.Measurement): the dx measure of the problem
+        """
         F = 0
         v = f.TestFunction(V)
-        self.post_processing_solution = f.Function(V)
-        F += -self.post_processing_solution * v * dx
+        c = f.TrialFunction(V)
+
+        F += -c * v * dx
         for mat in materials.materials:
             if mat.solubility_law == "sievert":
+                # for sievert materials c = theta * S
                 F += self.solution * self.S * v * dx(mat.id)
             elif mat.solubility_law == "henry":
+                # for henry materials c = theta**2 * S
                 F += self.solution**2 * self.S * v * dx(mat.id)
+
         self.form_post_processing = F
+        self.post_processing_solution = f.Function(V)