updated physics docs

docs/general/equilibria.md

@@ -8,7 +8,7 @@ last_modified_at: 2023-04-13
 ---
 
 For a complete list of implemented equilibria we refer to the [Legolas docs](../../ford/lists/modules.html), more specifically
-all descendants of the `mod_equilibrium` module.
+all descendants of the `mod_equilibrium` module. These all follow the naming convention `smod_equil_*` ("submodule equilibrium").
 You can take a look at the documentation of every submodule for information on default parameters, which ones you can vary and where
 they comes from. Below we give a small overview of which equilibria are implemented together with the default included
 physical effects, and links to the source docs.

docs/general/own_setup.md

@@ -5,8 +5,7 @@ classes: wide
 sidebar:
   nav: "leftcontents"
 toc: true
-toc_icon: "tools"
-toc_label: "User submodule configuration"
+toc_icon: "chevron-circle-down"
 last_modified_at: 2023-04-13
 ---
 
@@ -79,7 +78,7 @@ end module procedure user_defined_eq
 A third option is modifying the spacing function, which may be more convenient than defining a custom grid. An example is given below, where a Gaussian spacing function is used, resulting in grid accumulation near the peak.
 {% capture note %}
 <i class="fa fa-lightbulb" aria-hidden="true"></i>
-**Note:** the spacing function governs `dx`, meaning that smaller function values imply a smaller grid spacing, and hence more points. Also, when using this method the number of
+The spacing function governs `dx`, meaning that smaller function values imply a smaller grid spacing, and hence more points. Also, when using this method the number of
 gridpoints will vary to accomodate the specified spacing function.
 {% endcapture %}
 <div class="notice--info">
@@ -150,11 +149,11 @@ end procedure user_defined_eq
 **Disclaimer**: due to the finite element representation used by Legolas the equilibrium profiles (or their derivatives) do not have to be
 continuous. Do note though that smooth, continuous profiles give the best results, and you should try to achieve this as much as possible. Instead of a discontinuity
 you can try to connect different regions with a smooth sine or hyperbolic tangent profile. This is not at all required and Legolas will run just fine with strong jumps or gradients,
-but know that in that case you should double check if the results make sense and ensure enough gridpoints to resolve everything.
+but know that this may lead to artifacts in the eigenfunctions and/or spectrum.
 {: .notice--danger}
 
 Setting the actual equilibrium configuration is done by calling the corresponding type-bound procedure of the `background` object, and provide a function.
-All background functions take either no arguments, or a single `real(dp)` scalar argument indicating position. The example below sets a homogeneous density profile and a position-dependent temperature and $v_y$ profile.
+All background functions take either no arguments, or a single `real(dp)` scalar argument indicating position. The example below sets a homogeneous $\rho_0$ profile and a position-dependent $T_0$ and $v_y$ profile.
 For a full overview of the available background functions see [here](../../ford/type/background_t.html), or take a look at all the various
 pre-implemented equilibria [here](../../general/equilibria).
 ```fortran
@@ -174,8 +173,9 @@ contains
 
     ! Note: functions that are not set are automatically set to zero.
     call background%set_density_funcs(rho0_func=rho0)
-    call background%set_velocity_2_funcs(v02_func=v02, dv02_func=dv02, ddv02_func=ddv02)
+    ! dT0 denotes derivative with respect to position
     call background%set_temperature_funcs(T0_func=T0, dT0_func=dT0)
+    call background%set_velocity_2_funcs(v02_func=v02, dv02_func=dv02, ddv02_func=ddv02)
   end procedure user_defined_eq
 
   real(dp) function rho0()
@@ -210,7 +210,7 @@ Note that you can specify variables in module scope (like `alpha` in the example
 
 <i class="fas fa-lightbulb" aria-hidden="true"></i>
 **Tip:** to make sure your implementation is in order you can first run without solving the eigenvalue problem by setting `solver = "none"` in the `solvelist`.
-You can then plot the equilibrium profiles using Pylbo
+You can then plot the equilibrium profiles using Pylbo.
 {: .notice--info}
 
 
@@ -224,15 +224,15 @@ Below is an overview of the options for resistivity:
 ```fortran
 ! activate temperature-dependent resistivity (Spitzer)
 call settings%physics%enable_resistivity()
-
-! set resistivity to a constant value
+! OR set resistivity to a constant value
 call settings%physics%enable_resistivity(fixed_resistivity_value=0.0001_dp)
+
 ! additionally you can specify a dropoff profile
 settings%physics%resistivity%use_dropoff = .true.
 settings%physics%dropoff_edge_dist = 0.05_dp
 settings%physics%dropoff_width = 0.2_dp
 ```
-For a user-defined resistivity implementation you can repoint the corresponding function in the `physics` object, see [here](../../ford/type/physics_t.html) for more info.
+For a user-defined resistivity implementation you can repoint the corresponding function in the `physics` object, see [here](../../ford/type/physics_t.html#boundprocedure-set_resistivity_funcs) for more info.
 This is exactly the same as for the background functions, see an example below.
 
 ```fortran
@@ -277,22 +277,22 @@ In the absence of a magnetic field (hydrodynamics) perpendicular conduction uses
 call settings%physics%enable_parallel_conduction()
 ! activate perpendicular Spitzer conductivity
 call settings%physics%enable_perpendicular_conduction()
-
-! set conduction to a constant value
+! OR set conduction to a constant value
 call settings%physics%enable_parallel_conduction(fixed_tc_para_value=0.0001_dp)
 call settings%physics%enable_perpendicular_conduction(fixed_tc_perp_value=0.0001_dp)
 ```
-Both the parallel and perpendicular resistivity functions can be user-defined as well.
+Both the parallel and perpendicular resistivity functions can be user-defined as well using [these](../../ford/type/physics_t.html#boundprocedure-set_parallel_conduction_funcs) type-bound procedures.
 
 ### Radiative cooling and heating
 To have an equilibrium in thermal balance you should enable heating and leave the `force_thermal_balance` option to `.true.`.
 This will automatically set the heating in such a way that the thermal balance equation is satisfied.
-To run the code outside of thermal balance you can set this flag to false, and either disable heating alltogether or specify a user-defined heating function.
+To run the code outside of thermal balance you can set this flag to false.
 ```fortran
 call settings%physics%enable_cooling(cooling_curve="jc_corona")
 ! to achieve thermal balance, enable heating
 call settings%physics%enable_heating(force_thermal_balance=.true.)
 ```
+Both cooling and heating can be user-specified using [these](../../ford/type/physics_t.html#boundprocedure-set_cooling_funcs) type-bound procedures.
 
 ### Flow
 ```fortran
@@ -323,6 +323,7 @@ end submodule smod_user_defined
 ```fortran
 call settings%physics%enable_viscosity(viscosity_value=0.0001_dp, viscous_heating=.false.)
 ```
+Viscosity uses a constant value and does not accept a user-defined function.
 
 ### Hall MHD
 ```fortran
@@ -334,3 +335,4 @@ settings%physics%hall%use_inertia_dropoff
 settings%physics%dropoff_edge_dist = 0.05_dp
 settings%physics%dropoff_width = 0.2_dp
 ```
+The Hall treatment is fixed and does not accept user-defined functions.

docs/physics/equations.md

@@ -62,24 +62,27 @@ Resistivity can either be constant over the entire grid, or a more general tempe
 by the Spitzer resistivity:
 
 $$
-\eta = \frac{4\sqrt{2\pi}}{3}\frac{Z_\text{ion}e^2\sqrt{m_e}\ln(\lambda)}{(4\pi\epsilon_0)^2(k_B T)^{3/2}},
+\eta(T) = \frac{4\sqrt{2\pi}}{3}\frac{Z_\text{ion}e^2\sqrt{m_e}\ln(\lambda)}{(4\pi\epsilon_0)^2(k_B T)^{3/2}},
 $$
 
-with the Coulomb logarithm $\ln(\lambda) \approx 22$.
+with the Coulomb logarithm $\ln(\lambda) \approx 22$. The function $\eta(T)$ can be user-specified, in which case its temperature derivative should be provided as well.
+If the resistivity profile explicitly depends on position as well, providing the derivative of $\eta(x, T)$ with respect to $x$ is an additional requirement.
 
-### Optically thin radiative losses
+### Heating and optically thin radiative losses
 Radiative cooling is governed by the heat-loss function, specified as the difference between energy gains and energy losses
 
 $$
-\HL = \rho\Lambda(T) - \HH,
+\HL = \rho\Lambda(T) - \HH(\rho, T),
 $$
 
 The function $\Lambda(T)$ here is called the _cooling curve_, which is a tabulated set of values resulting from detailed molecular calculations.
 Legolas has multiple cooling curves from existing literature implemented which are interpolated at high resolution and sampled on the given grid.
 We have also included analytical, piecewise prescriptions (e.g. the `rosner` curve). See the [physicslist](../../general/parameter_file/#physicslist) for an overview of the different options.
+The function $\Lambda(T)$ can be user-specified, in which case its temperature derivative should be provided as well.
 
-The function $\HH$ is called the _heating term_, and it is currently not (yet) possible to specify this. We assume that this term only depends on the equilibrium background, meaning that it is
-constant in time but possibly varying in space and as such that it balances out the cooling contribution to ensure thermal equilibrium.
+The function $\HH(\rho, T)$ specifies the heating function. If thermal balance is forced we assume that this term only depends on the equilibrium background, meaning that it is
+constant in time but possibly varying in space and as such that it balances out the cooling contribution to ensure thermal equilibrium. $\HH(\rho, T)$ can be user-specified, in which case its density and
+temperature derivatives should be provided as well.
 
 ### Thermal conduction
 The thermal conduction prescription relies on a tensor representation to model the anisotropy in MHD and is given by
@@ -98,6 +101,25 @@ $$
 $$
 
 Alternatively, both of these can be independently set to constant values.
+In the absence of a magnetic field we set $\kappa_\bot = \kappa_\parallel$, such that the tensor representation reduces to $\kappabf = \kappa_\bot\boldsymbol{I}$.
+
+When overriding parallel thermal conduction a function $\kappa_\parallel(T)$ can be given, as well as its temperature derivative. The derivative with respect to position (denoted with $'$) is automatically calculated as $\kappa_\parallel' = \frac{\partial \kappa_\parallel}{\partial T}T_0'$.
+
+For perpendicular thermal conduction a function $\kappa_\bot(\rho, T, B)$ can be given, as well as its derivatives
+
+$$
+\frac{\partial \kappa_\bot}{\partial \rho}, \qquad
+\frac{\partial \kappa_\bot}{\partial T}, \qquad
+\frac{\partial \kappa_\bot}{\partial B^2}.
+$$
+
+The derivative with respect to position is automatically calculated as
+
+$$ \kappa_\bot' =
+\frac{\partial \kappa_\bot}{\partial T}T_0'
++ \frac{\partial \kappa_\bot}{\partial \rho}\rho_0'
++ \frac{\partial \kappa_\bot}{\partial B^2}2B_0B_0'.
+$$
 
 ### Viscosity
 The viscosity addition to the momentum equation is usually given in a full tensor form, but in Legolas this is simplified to good approximation to

docs/physics/units.md

@@ -4,72 +4,69 @@ layout: single
 classes: wide
 sidebar:
   nav: "leftcontents"
-last_modified_at: 2021-07-27
+last_modified_at: 2023-04-23
 ---
 
 All equations in Legolas are in dimensionless form, as is common practice when dealing with (M)HD.
-
-As usual we have **three** degrees of freedom: the unit magnetic field $B_u$ and unit length $L_u$ should be set,
-together with _either_ the unit density $\rho_u$ **OR** unit temperature $T_u$.
+As usual we have **three** degrees of freedom.
 
 ## Mean molecular weight
 Unit normalisations depend on the molecular weight $\bar{\mu}$, and in Legolas we usually distinguish between two cases:
 
-- **Pure proton plasma**: $\bar{\mu} = 1$, this is the default case.
+- **Electron-proton plasma**: $\bar{\mu} = 0.5$, this is the default case.
 
     $$
-    \bar{\mu} = \dfrac{m_i n_i}{n_i} = m_i \rightarrow \bar{\mu} = 1
+    \bar{\mu} = \dfrac{m_e n_e + m_i n_i}{n_e + n_i} \simeq \dfrac{m_i n_i}{n_e + n_i}
+              = \dfrac{1}{2}m_i \rightarrow \bar{\mu} = \dfrac{1}{2}
     $$
 
-- **Electron-proton plasma**: $\bar{\mu} = 0.5$, in this case the molecular weight should be explicitly set.
+- **Pure proton plasma**: $\bar{\mu} = 1$, in this case the molecular weight should be explicitly set.
 
     $$
-    \bar{\mu} = \dfrac{m_e n_e + m_i n_i}{n_e + n_i} \simeq \dfrac{m_i n_i}{n_e + n_i}
-              = \dfrac{1}{2}m_i \rightarrow \bar{\mu} = \dfrac{1}{2}
+    \bar{\mu} = \dfrac{m_i n_i}{n_i} = m_i \rightarrow \bar{\mu} = 1
     $$
 
 ## Normalisations
-We have two options: either the unit density is set or the unit temperature is set. In both cases a unit magneticfield
-and a unit length should be set as well, so the unit pressure is given by
+Legolas has three options to specify units, all in cgs. In what follows $m_p$ denotes the proton mass,
+$k_B$ the Boltzmann constant, and $\mu_0 = 4\pi$ the magnetic constant.
 
-$$
-p_u = \dfrac{B_u^2}{\mu_0},
-$$
+1. Reference unit density, unit magnetic field and unit length $(\rho_u, B_u, L_u)$, then
 
-taking a plasma-beta reference value of 2. The magnetic constant is given by $\mu_0$ ($= 4\pi$ in cgs).
+   $$
+   p_u = \frac{B_u^2}{\mu_0}, \quad
+   T_u = \frac{\bar{\mu} p_u m_p}{k_B \rho_u}, \quad
+   n_u = \frac{\rho_u}{m_p}, \quad
+   v_u = \frac{B_u}{\sqrt{\mu_0 \rho_u}}.
+   $$
 
-We then have the following options:
+2. Reference unit temperature, unit magnetic field and unit length $(T_u, B_u, L_u)$, then
 
-- If $\rho_u$ chosen along with $B_u$ and $L_u$:
+   $$
+   p_u = \frac{B_u^2}{\mu_0}, \quad
+   \rho_u = \frac{\bar{\mu} p_u m_p}{k_B T_u}, \quad
+   n_u = \frac{\rho_u}{m_p}, \quad
+   v_u = \frac{B_u}{\sqrt{\mu_0 \rho_u}}.
+   $$
 
-    $$
-    T_u = \dfrac{\bar{\mu}m_p p_u}{k_B \rho_u},
-    $$
-
-- If $T_u$ chosen along with $B_u$ and $L_u$:
-
-    $$
-    \rho_u = \dfrac{\bar{\mu}m_p p_u}{k_B T_u},
-    $$
+3. Reference unit numberdensity, unit temperature and unit length $(n_u, T_u, L_u)$, then
 
-All other normalisations follow from these, where the Alfvén speed is assumed as reference velocity:
+   $$
+   p_u = \bar{\mu} n_u k_B T_u, \quad
+   \rho_u = m_p n_u, \quad
+   v_u = \sqrt{\frac{p_u}{\rho_u}}, \quad
+   B_u = \sqrt{\mu_0 p_u}.
+   $$
 
-$$
-\begin{gather}
-    m_u = \rho_u L_u, \\
-    n_u = \dfrac{\rho_u}{m_p}, \\
-    v_u = \dfrac{B_u}{\sqrt{\mu_0 \rho_u}}, \\
-    t_u = \dfrac{L_u}{v_u}, \\
-    \Lambda_u = \dfrac{p}{t_u n_u^2}, \\
-    \kappa_u = \dfrac{\rho_u L_u v_u^3}{T_u},
-\end{gather}
-$$
+All other normalisations follow from those above and are given by
+- unit mass: $M_u = \rho_u L_u^3$
+- unit time: $t_u = \dfrac{L_u}{v_u}$
+- unit resistivity: $\eta_u = \dfrac{L_u^2}{t_u}$
+- unit cooling curve: $\Lambda_u = \dfrac{p_u}{t_u n_u^2}$
+- unit conduction: $\kappa_u = \dfrac{\rho_u L_u v_u^3}{T_u}$
 
-which, from top to bottom, read as mass unit, numberdensity unit, velocity unit, time unit, cooling curve unit and
-conduction unit. The proton mass is given by $m_p$ and $k_B$ denotes the Boltzmann constant.
 
 <i class="fas fa-lightbulb" aria-hidden="true"></i>
-**Note:** The unit normalisations are only relevant when radiative cooling, thermal conduction or temperature-dependent
-resistivity is included. We always set base values though (as one should), which are given in cgs units by
-$B_u = 10$ G, $L_u = 10^9$ cm and $T_u = 10^6$ K. If no normalisations are specified these are used by default.
-{: .notice--info}
+**Note:** the unit normalisations are only relevant when radiative cooling, thermal conduction or temperature-dependent resistivity is included.
+We always set base values though (as one should), which are set using option 2. with default values
+$B_u = 10$ G, $L_u = 10^9$ cm and $T_u = 10^6$ K.
+{: .notice--info}

docs/pylbo/legolas_interface.md

@@ -5,7 +5,7 @@ classes: wide
 sidebar:
   nav: "leftcontents"
 toc: true
-last_modified_at: 2021-07-27
+last_modified_at: 2021-04-13
 ---
 {% capture tip %}
 <i class="fas fa-lightbulb" aria-hidden="true"></i>
@@ -142,8 +142,7 @@ config = {
     },
     "equilibrium_type": "resistive_homo",
     "resistivity": True,
-    "use_fixed_resistivity": True,
-    "fixed_eta_value": np.linspace(0.001, 0.01, 50),
+    "fixed_resistivity_value": np.linspace(0.001, 0.01, 50),
     "logging_level": 0,
     "show_results": False,
     "write_eigenfunctions": True,
@@ -170,8 +169,8 @@ you can supply a full path instead, either as a string (e.g.`output_dir="/users/
 ## Running Legolas with Pylbo
 The parfiles generated in the above examples can be passed on to Pylbo, which in turn will pass those on to Legolas.
 
-**Remark**: Note that Pylbo will always link to the executable in the Legolas home directory by default. It is therefore good
-practice to explicitly specify the Legolas executable when calling the runner.
+<i class="fa fa-exclamation-triangle" aria-hidden="true"></i>
+Note that Pylbo will always link to the executable in the Legolas home directory by default. It is therefore good practice to explicitly specify the Legolas executable when calling the runner.
 {: .notice--danger}
 
 ### Single run