Finally! Dividing by the volume of the node gave the right value

src/modules/common/moduleRefParam.f90

@@ -5,5 +5,35 @@ MODULE moduleRefParam
   !Reference parameters for non-dimensional problem
   REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
 
+  contains
+    subroutine calculateDerivedRef()
+      use moduleConstParam, only: pi, qe, eps_0, kb
+      implicit none
+      
+      !Derived parameters
+      L_ref     = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
+      v_ref     = DSQRT(kb*T_ref/m_ref)          !reference velocity
+      ti_ref    = L_ref/v_ref                    !reference time
+      Vol_ref   = L_ref**3                       !reference volume
+      EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
+      Volt_ref  = EF_ref*L_ref                   !reference voltage
+      B_ref     = m_ref / (ti_ref * qe)          !reference magnetic field
+
+      if (sigmaVrel_ref == 0.0d0) then
+        ! No reference cross section provided, check radius
+        if (r_ref == 0.0d0) then
+          ! No reference radius provided, use reference length
+          sigmaVrel_ref = L_ref**2 * v_ref
+
+        else
+          ! Use the reference radius
+          sigmaVrel_ref = pi*(r_ref+r_ref)**2*v_ref
+
+        end if
+
+      end if
+
+    end subroutine calculateDerivedRef
+
 END MODULE moduleRefParam
 

src/modules/init/moduleInput.f90

@@ -111,7 +111,6 @@ MODULE moduleInput
     !Reads the reference parameters
     SUBROUTINE readReference(config)
       USE moduleRefParam
-      USE moduleConstParam
       USE moduleErrors
       USE json_module
       IMPLICIT NONE
@@ -131,30 +130,14 @@ MODULE moduleInput
       CALL config%get(object // '.temperature', T_ref, found)
       IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
 
-      !Derived parameters
-      L_ref     = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
-      v_ref     = DSQRT(kb*T_ref/m_ref)          !reference velocity
-      ti_ref    = L_ref/v_ref                    !reference time
-      Vol_ref   = L_ref**3                       !reference volume
-      EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
-      Volt_ref  = EF_ref*L_ref                   !reference voltage
-      B_ref     = m_ref / (ti_ref * qe)          !reference magnetic field
+      ! initialize values for cross sections
+      sigmaVrel_ref = 0.0d0
+      r_ref = 0.0d0
 
-      !If a reference cross section is given, it is used
       CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
-
-      !If not, the reference radius is searched
-      IF (.NOT. found) THEN
       CALL config%get(object // '.radius', r_ref, found)
-        IF (found) THEN
-          sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity
 
-        ELSE
-          sigmaVrel_ref = L_ref**2 * v_ref
-
-        END IF
-
-      END IF
+      call calculateDerivedRef()
 
     END SUBROUTINE readReference
 

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -163,7 +163,7 @@ submodule(moduleMesh) boundaryEM
 
       do n = 1, self%nNodes
         associate(node => self%nodes(n)%obj)
-          vectorF(node%n) = vectorF(node%n) + self%electricField
+          vectorF(node%n) = vectorF(node%n) - self%electricField/node%v
 
         end associate