Fix signs in advection-diffusion equations and update comments accordingly

Author: nikoscham

examples/generalFormPDEScript/advectionDiffusion1D/README.md

@@ -5,7 +5,7 @@
 This example demonstrates solving a one-dimensional advection-diffusion problem with a Gaussian source term using the FEAScript library. The problem models the transport equation:
 
 $$
-\frac{d^2u}{dx^2} + 10 \frac{du}{dx} = -10 \cdot e^{-200 \cdot (x - 0.5)^2}
+\frac{d^2u}{dx^2} - 10 \frac{du}{dx} = 10 \cdot e^{-200 \cdot (x - 0.5)^2}
 $$
 
 This can be written in the general form as:
@@ -17,9 +17,9 @@ $$
 where:
 
 - $A(x) = 1$ (diffusion coefficient)
-- $B(x) = 10$ (advection coefficient)
+- $B(x) = -10$ (advection coefficient)
 - $C(x) = 0$ (no reaction term)
-- $D(x) = -10 \cdot e^{-200 \cdot (x - 0.5)^2}$ (Gaussian source term centered at x = 0.5)
+- $D(x) = 10 \cdot e^{-200 \cdot (x - 0.5)^2}$ (Gaussian source term centered at x = 0.5)
 
 ### Instructions
 

examples/generalFormPDEScript/advectionDiffusion1D/advectionDiffusion1D.js

@@ -23,11 +23,11 @@ const model = new FEAScriptModel();
 // Set solver configuration with coefficient functions
 model.setSolverConfig("generalFormPDEScript", {
   coefficientFunctions: {
-    // Equation d²u/dx² + 10 du/dx = -10 * exp(-200 * (x - 0.5)²)
+    // Equation d²u/dx² - 10 du/dx = 10 * exp(-200 * (x - 0.5)²)
     A: (x) => 1, // Diffusion coefficient
-    B: (x) => 10, // Advection coefficient
+    B: (x) => -10, // Advection coefficient
     C: (x) => 0, // Reaction coefficient
-    D: (x) => -10 * Math.exp(-200 * Math.pow(x - 0.5, 2)), // Source term
+    D: (x) => 10 * Math.exp(-200 * Math.pow(x - 0.5, 2)), // Source term
   },
 });
 

src/solvers/generalFormPDEScript.js

@@ -95,7 +95,7 @@ export function assembleGeneralFormPDEMat(meshData, boundaryConditions, coeffici
           const globalNodeIndex1 = localToGlobalMap[localNodeIndex1];
 
           // Source term contribution to residual vector
-          residualVector[globalNodeIndex1] +=
+          residualVector[globalNodeIndex1] -=
             gaussWeights[gaussPointIndex] * detJacobian * d * basisFunction[localNodeIndex1];
 
           for (let localNodeIndex2 = 0; localNodeIndex2 < numNodes; localNodeIndex2++) {
@@ -110,15 +110,15 @@ export function assembleGeneralFormPDEMat(meshData, boundaryConditions, coeffici
               basisFunctionDerivX[localNodeIndex2];
 
             // Advection term
-            jacobianMatrix[globalNodeIndex1][globalNodeIndex2] +=
+            jacobianMatrix[globalNodeIndex1][globalNodeIndex2] -=
               gaussWeights[gaussPointIndex] *
               detJacobian *
               b *
               basisFunctionDerivX[localNodeIndex2] *
               basisFunction[localNodeIndex1];
 
             // Reaction term
-            jacobianMatrix[globalNodeIndex1][globalNodeIndex2] -=
+            jacobianMatrix[globalNodeIndex1][globalNodeIndex2] +=
               gaussWeights[gaussPointIndex] *
               detJacobian *
               c *
@@ -219,7 +219,7 @@ export function assembleGeneralFormPDEFront({
       // Computation of local Jacobian matrix and residual vector
       for (let localNodeIndex1 = 0; localNodeIndex1 < numNodes; localNodeIndex1++) {
         // Source term contribution to local residual vector
-        localResidualVector[localNodeIndex1] +=
+        localResidualVector[localNodeIndex1] -=
           gaussWeights[gaussPointIndex] * detJacobian * d * basisFunction[localNodeIndex1];
 
         for (let localNodeIndex2 = 0; localNodeIndex2 < numNodes; localNodeIndex2++) {
@@ -232,15 +232,15 @@ export function assembleGeneralFormPDEFront({
             basisFunctionDerivX[localNodeIndex2];
 
           // Advection term
-          localJacobianMatrix[localNodeIndex1][localNodeIndex2] +=
+          localJacobianMatrix[localNodeIndex1][localNodeIndex2] -=
             gaussWeights[gaussPointIndex] *
             detJacobian *
             b *
             basisFunctionDerivX[localNodeIndex2] *
             basisFunction[localNodeIndex1];
 
           // Reaction term
-          localJacobianMatrix[localNodeIndex1][localNodeIndex2] -=
+          localJacobianMatrix[localNodeIndex1][localNodeIndex2] +=
             gaussWeights[gaussPointIndex] *
             detJacobian *
             c *