parallel_ohms_law: Added a new component. by malamast · Pull Request #432 · boutproject/hermes-3

malamast · 2025-12-02T22:42:14Z

Uses parallel Ohm's law to calculate the parallel current
Then uses this parallel current to calculate the parallel electron velocity.
It can potentially help with numerical stability in 2D asymmetric simulations (Tokamak plasma edge sims) when drifts are included.
This component requires collisions to be calculated first, so it should be called after collisions component.
When parallel_ohms_law is used then there is no need to evolve_momentum for electrons

/// Caution! This component is not carfully evaluated for multispecies cases with multiple ions.

/// Use parallel Ohm's law to calculate the parallel current.
/// This law basically comes from the electron parallel momentum equation
/// by neglecting inertial (small by m_e / m_i) and viscous terms
/// (see Phys. Plasmas 10, 4744–4757 (2003) by Simakov & Catto).
/// Use the parallel current to calculate the electron parallel velocity
///
/// Jpar = - 1/eta * ∇phi + 1/(eta e n_e) * ∇P_e + (0.71 k_B)/(eta e) * ∇T_e
///
/// where phi is the electric potential and eta is the resistivity.
///
/// After normalization the equation takes the form:
///
/// Jpar = - ∇phi / eta + ∇P_e / (eta n_e) + 0.71 ∇T_e / eta
///
/// where the normalization parameter for the resistivity eta is:
///
/// eta_norm = Tnorm / (rho_s0 * Nnorm * Cs0 * qe)
///
/// Then uses this parallel current to calculate the parallel electron velocity.
///
/// Note: This needs to be called after collisions and other
/// components which impose forces on electrons_

codecov · 2025-12-02T22:49:15Z

Codecov Report

❌ Patch coverage is 0% with 170 lines in your changes missing coverage. Please review.
✅ Project coverage is 22.03%. Comparing base (0e601f7) to head (c8a1997).

Files with missing lines	Patch %	Lines
src/parallel_ohms_law.cxx	0.00%	167 Missing ⚠️
include/parallel_ohms_law.hxx	0.00%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #432      +/-   ##
==========================================
- Coverage   22.50%   22.03%   -0.47%     
==========================================
  Files          87       89       +2     
  Lines        8071     8241     +170     
  Branches     1146     1167      +21     
==========================================
  Hits         1816     1816              
- Misses       6085     6255     +170     
  Partials      170      170

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

mikekryjak · 2025-12-03T09:08:55Z

Hi @malamast, nice work. Can you tell me a bit more about how mean field drifts with this approach compare to what SOLPS, UEDGE etc are doing? It looks more or less like what I understand is in SOLPS, but I'm not a drifts expert and never looked at the equations in detail.

ZedThree

Some little bits to tidy up, but otherwise looks good. I'll leave the physics review up to @mikekryjak / @bendudson !

ZedThree · 2025-12-03T10:00:26Z

+
+
+namespace {
+BoutReal clip(BoutReal value, BoutReal min, BoutReal max) {


This is unused and can be deleted

ZedThree · 2025-12-03T10:00:46Z

+  return value;
+}
+
+Ind3D indexAt(const Field3D& f, int x, int y, int z) {


This is in hermes_utils, so this implementation can be deleted

ZedThree · 2025-12-03T10:00:57Z

+///
+/// exp( 2*log(fc) - log(fm) )
+///
+BoutReal limitFree(BoutReal fm, BoutReal fc) {


Also unused and can be deleted here

ZedThree · 2025-12-03T10:02:33Z

+  }
+
+  eta = softFloor( resistivity_eta / eta_norm , resistivity_floor );
+  // eta = sqrt( SQ(resistivity_eta / eta_norm) + resistivity_floor*resistivity_floor ); // This mitigates discontinuities.


We can delete commented out code:

Suggested change

// eta = sqrt( SQ(resistivity_eta / eta_norm) + resistivity_floor*resistivity_floor ); // This mitigates discontinuities.

ZedThree · 2025-12-03T10:02:53Z

+      // if (!current.isAllocated()) {
+        // Not yet allocated -> Set to the value
+        // This avoids having to set to zero initially and add the first time
+        // current = charge * N * V;
+      // } else {
+      current += charge * N * V;
+      // }


Suggested change

// if (!current.isAllocated()) {

// Not yet allocated -> Set to the value

// This avoids having to set to zero initially and add the first time

// current = charge * N * V;

// } else {

current += charge * N * V;

// }

current += charge * N * V;

ZedThree · 2025-12-03T10:03:37Z

+  // set(electrons["velocity"], Ve); // # This is causing a problem with option.hasAttribute("final-domain") check in set function.
+  electrons["velocity"].force(Ve);


I'm not quite sure what this comment means -- is there an issue here that needs to be resolved?

The problem is probably that:

This component needs the collision frequency, so needs to run after collisions

collisions uses the species velocities so marks electron velocity as final

This component needs to set the velocity. This isn't consistent with the friction calculation in collisions so an exception is thrown. This is the mechanism working as intended, flagging a circular dependency.

I think this should be resolved by #399, because that separates the collision frequency calculation from the friction calculation. This parallel_ohms_law component could run between braginskii_collisions and braginskii_friction.

1) Use parallel Ohm's law to calculate the parallel current 2) then uses this parallel current to calculate the parallel electron velocity. -This compoment requries cossisions to be calculated first so it should be called after collisions component. -When parallel_ohms_law is used then there is no need to evolve_momentum for electrons

…sed methods.

malamast · 2025-12-23T00:38:57Z

Hi @malamast, nice work. Can you tell me a bit more about how mean field drifts with this approach compare to what SOLPS, UEDGE etc are doing? It looks more or less like what I understand is in SOLPS, but I'm not a drifts expert and never looked at the equations in detail.

In SOLPS and UEDGE, the parallel Ohm’s law is used to compute the electron parallel velocity, rather than solving a full electron momentum equation. I implemented this approach with the expectation that it would improve the numerical stability of my 2D transport simulations when drift terms are included.

In practice, I enable this through the following component setup:
components = (d+, d, e,
collisions, relax_potential,
sheath_boundary_simple, parallel_ohms_law,
recycling, reactions,
ion_viscosity,
diamagnetic_drift, polarisation_drift,
thermal_force, sound_speed,
fixed_fraction_hutchinson_carbon)

For the electron species, the configuration is:
[e]
type = quasineutral, evolve_pressure, anomalous_diffusion

Based on some runs with this formulation, the resulting electron parallel velocity Ve was lower than when using evolve_momentum but the trends were the same (except some initial transient phases with the evolve_momentum).

However, at present, I do not have a reliable validation strategy for this component because I couldn't make it work with the drifts. I suggest not merging it for now.

…d before the sheath_boundary component. using IS_SET_NOBOUNDARY and GET_NOBOUNDARY macros so that it can set Ve and NVe in the electron state without forcing it.

malamast · 2026-01-10T01:51:33Z

This is in better shape now. It needs to be called after braginskii_collisions and after vorticity or relax_potential.
Note that when we use parallel_ohms_law then we do not use evolve_momentum for electrons.

  components = (d+, e, d,
                braginskii_collisions, relax_potential,
                parallel_ohms_law, sheath_boundary_simple,
                braginskii_friction, braginskii_heat_exchange,
                sound_speed, reactions,
                braginskii_ion_viscosity, braginskii_electron_viscosity,
                diamagnetic_drift,
                braginskii_thermal_force, braginskii_conduction, recycling)

  [e]
  type = quasineutral, evolve_pressure, anomalous_diffusion

I'm preparing some comparisons with the evolve_momentum case.

malamast requested review from ZedThree, bendudson and mikekryjak December 2, 2025 22:42

malamast assigned bendudson Dec 2, 2025

malamast added enhancement New feature or request help wanted Extra attention is needed labels Dec 2, 2025

ZedThree requested changes Dec 3, 2025

View reviewed changes

malamast added 2 commits December 22, 2025 16:21

parallel_ohms_law: Tidied up the code by removing unnecessary and unu…

d421c5a

…sed methods.

malamast force-pushed the parallel_ohms_law branch from c8a1997 to d421c5a Compare December 23, 2025 00:24

malamast added 2 commits January 9, 2026 14:56

Removed some old renamimng conflicts related to electron_viscosity.

cd35174

parallel_ohms_law: Restractured the component so that it can be calle…

1e3bbe0

…d before the sheath_boundary component. using IS_SET_NOBOUNDARY and GET_NOBOUNDARY macros so that it can set Ve and NVe in the electron state without forcing it.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

parallel_ohms_law: Added a new component.#432

parallel_ohms_law: Added a new component.#432
malamast wants to merge 4 commits into
masterfrom
parallel_ohms_law

malamast commented Dec 2, 2025

Uh oh!

codecov Bot commented Dec 2, 2025

Uh oh!

mikekryjak commented Dec 3, 2025 •

edited

Loading

Uh oh!

ZedThree left a comment

Uh oh!

ZedThree Dec 3, 2025

Uh oh!

ZedThree Dec 3, 2025

Uh oh!

ZedThree Dec 3, 2025

Uh oh!

ZedThree Dec 3, 2025

Uh oh!

ZedThree Dec 3, 2025

Uh oh!

ZedThree Dec 3, 2025

Uh oh!

bendudson Dec 9, 2025

Uh oh!

malamast commented Dec 23, 2025

Uh oh!

malamast commented Jan 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants



		namespace {
		BoutReal clip(BoutReal value, BoutReal min, BoutReal max) {

		// set(electrons["velocity"], Ve); // # This is causing a problem with option.hasAttribute("final-domain") check in set function.
		electrons["velocity"].force(Ve);

Conversation

malamast commented Dec 2, 2025

Uh oh!

codecov Bot commented Dec 2, 2025

Codecov Report

Uh oh!

mikekryjak commented Dec 3, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ZedThree left a comment

Choose a reason for hiding this comment

Uh oh!

ZedThree Dec 3, 2025

Choose a reason for hiding this comment

Uh oh!

ZedThree Dec 3, 2025

Choose a reason for hiding this comment

Uh oh!

ZedThree Dec 3, 2025

Choose a reason for hiding this comment

Uh oh!

ZedThree Dec 3, 2025

Choose a reason for hiding this comment

Uh oh!

ZedThree Dec 3, 2025

Choose a reason for hiding this comment

Uh oh!

ZedThree Dec 3, 2025

Choose a reason for hiding this comment

Uh oh!

bendudson Dec 9, 2025

Choose a reason for hiding this comment

Uh oh!

malamast commented Dec 23, 2025

Uh oh!

malamast commented Jan 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

mikekryjak commented Dec 3, 2025 •

edited

Loading