[PL/CT] compute A and rhs.

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -278,13 +278,10 @@ void ComponentTransportProcess::solveReactionEquation(
         matrix_specification, matrix_id);
     auto& b =
         NumLib::GlobalVectorProvider::provider.getVector(matrix_specification);
-    auto& rhs =
-        NumLib::GlobalVectorProvider::provider.getVector(matrix_specification);
 
     M.setZero();
     K.setZero();
     b.setZero();
-    rhs.setZero();
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &ComponentTransportLocalAssemblerInterface::assembleReactionEquation,
@@ -296,17 +293,30 @@ void ComponentTransportProcess::solveReactionEquation(
     MathLib::finalizeMatrixAssembly(K);
     MathLib::finalizeVectorAssembly(b);
 
-    MathLib::LinAlg::scale(M, 1.0 / dt);
-    MathLib::LinAlg::matMultAdd(M, *x[process_id], b, rhs);
+    auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(
+        matrix_specification, matrix_id);
+    auto& rhs =
+        NumLib::GlobalVectorProvider::provider.getVector(matrix_specification);
+
+    A.setZero();
+    rhs.setZero();
+
+    // compute A
+    MathLib::LinAlg::copy(M, A);
+    MathLib::LinAlg::aypx(A, 1.0 / dt, K);
+
+    // compute rhs
+    MathLib::LinAlg::matMult(M, *x[process_id], rhs);
+    MathLib::LinAlg::aypx(rhs, 1.0 / dt, b);
 
     using Tag = NumLib::NonlinearSolverTag;
     using EqSys = NumLib::NonlinearSystem<Tag::Picard>;
     auto& equation_system = static_cast<EqSys&>(ode_sys);
-    equation_system.applyKnownSolutionsPicard(M, rhs, *x[process_id]);
+    equation_system.applyKnownSolutionsPicard(A, rhs, *x[process_id]);
 
     auto& linear_solver =
         _process_data.chemical_solver_interface->linear_solver;
-    linear_solver.solve(M, rhs, *x[process_id]);
+    linear_solver.solve(A, rhs, *x[process_id]);
 
     NumLib::GlobalMatrixProvider::provider.releaseMatrix(M);
     NumLib::GlobalMatrixProvider::provider.releaseMatrix(K);