Polishing the docs

docs/src/man/analysis.md

@@ -96,7 +96,7 @@ Epx, Epy, Epz = get_pe_E(bd, species=0)
 
 The `get_pe_E` function calculates the divergence of the electron pressure tensor. For structured grids, it uses central finite differences. The species mass is automatically retrieved from simulation parameters (e.g., `ms0`, `mass0`), but can be overridden with the `mass` keyword argument. Note that for `PLANETARY` or `km` units, the resulting electric field is scaled to $\mu V/m$ for consistency with the convection and Hall terms.
 
-For structured grids, the current density ``\\mathbf{J} = \\nabla \\times \\mathbf{B}`` can also be calculated using finite differences with `get_current_density`, which returns a tuple of `(Jx, Jy, Jz)`.
+For structured grids, the current density $\mathbf{J} = \nabla \times \mathbf{B}$ can also be calculated using finite differences with `get_current_density`, which returns a tuple of `(Jx, Jy, Jz)`.
 
 ```julia
 Jx, Jy, Jz = get_current_density(bd)
@@ -119,12 +119,12 @@ Here is a full list of predefined derived quantities available via symbol indexi
 | :U                   | Velocity vector                  | Ux, Uy, Uz         |
 | :U0                  | Species 0 velocity vector        | UxS0, UyS0, UzS0   |
 | :U1                  | Species 1 velocity vector        | UxS1, UyS1, UzS1   |
-| :J                   | Current density vector (``\\mu A/m^2``)* | Bx, By, Bz         |
+| :J                   | Current density vector ($\mu A/m^2$)* | Bx, By, Bz         |
 | :b, :b2              | Magnetic field magnitude (sq)    | Bx, By, Bz         |
 | :e                   | Electric field magnitude         | Ex, Ey, Ez         |
 | :u                   | Velocity magnitude               | Ux, Uy, Uz         |
-| :j                   | Current density magnitude (``\\mu A/m^2``)* | Bx, By, Bz         |
-| :jx, :jy, :jz        | Current density components (``\\mu A/m^2``)* | Bx, By, Bz         |
+| :j                   | Current density magnitude ($\mu A/m^2$)* | Bx, By, Bz         |
+| :jx, :jy, :jz        | Current density components ($\mu A/m^2$)* | Bx, By, Bz         |
 | :anisotropy0, 1, ... | Pressure anisotropy              | B and P tensor     |
 
 \* Units are for PLANETARY data. For other systems, the units depend on the simulation normalization.
@@ -134,7 +134,17 @@ Here is a full list of predefined derived quantities available via symbol indexi
 To get the index of a certain quantity, e.g. electron number density
 
 ```julia
-ρe_= findfirst(x->x=="rhos0", bd.head.wname)
+id = findindex(bd, "rhos0")
+```
+
+#### Get grid range
+
+To extract the coordinates of the simulation grid:
+
+```julia
+ranges = get_range(bd)
+# For 2D, this is equivalent to:
+xrange, yrange = get_range(bd)
 ```
 
 #### Get variable range

docs/src/man/log.md

@@ -17,6 +17,16 @@ datapath = artifact"testdata" # where you can find multiple test data files
 
 These are also used in the standard test. These will be automatically downloaded from [batsrus_data](https://github.com/henry2004y/batsrus_data) if you run the package test locally.
 
+### BATL AMR Grid Structure
+
+The Block Adaptive Tree Library (BATL) used in Batsrus implements a block-based non-overlapping AMR scheme. Its key characteristics are:
+
+1. **Block-Based**: The computational domain is divided into blocks (or patches). Every block contains exactly the same number of grid cells (e.g., $8 \times 8 \times 8$ in 3D).
+2. **Octree Refinement**: When a region requires higher resolution, an entire block is refined by splitting it into $2^{nDim}$ child blocks ($8$ in 3D, $4$ in 2D). Each child block again has the same number of cells as the parent, meaning the spatial resolution is effectively doubled in each dimension.
+3. **Non-Overlapping**: Only the leaf nodes of the AMR tree (the blocks at the highest local level of refinement) are active and used for computations. Child blocks completely replace their parent blocks. Therefore, active grids at different levels never overlap, and active grids at the same level abut each other perfectly.
+
+In the `src/vtk.jl` code, the BATL AMR tree is managed via the `iTree_IA` array, which keeps track of each node's status (used, unused, refine, coarsen), level, local block index, spatial coordinates, and parent/child relationships.
+
 ## VTK AMR Grid Structure
 
 `vtkOverlappingAMR` implements a somewhat strict Berger-Collela AMR scheme:
@@ -66,4 +76,21 @@ If we can directly tell ParaView that the mesh we have is a dual-mesh, then the
 
 `AMRGaussianPulseSource`
 
-See [Multi-Resolution Rendering with Overlapping AMR](https://www.paraview.org/ParaView/index.php/Multi-Resolution_Rendering_with_Overlapping_AMR) for the implementation of C++ reader in VTK.
+See [Multi-Resolution Rendering with Overlapping AMR](https://www.paraview.org/ParaView/index.php/Multi-Resolution_Rendering_with_Overlapping_AMR) for the implementation of C++ reader in VTK.
+
+## Investigation: Converting BATL AMR to VTK AMR Formats
+
+Currently, `BatsrusWriteVTKExt` converts the BATL AMR grid into an unstructured VTK format (`vtkUnstructuredGrid`) by flattening the active blocks and writing them out as an array of `VTK_QUAD` (2D) or `VTK_HEXAHEDRON` (3D) cells. 
+
+### Limitations of the Current Unstructured Approach
+* **Loss of Hierarchy**: The multiresolution tree structure is discarded. VTK/ParaView cannot take advantage of the AMR levels for optimized rendering or dual-grid contouring.
+* **Overhead**: Unstructured grids explicitly store connectivity for every single cell. Since BATL blocks are Cartesian and structured internally, storing explicit cell connectivity is highly redundant and leads to unnecessarily large file sizes and increased memory consumption.
+
+### Better Approaches
+A more optimal approach is to map the BATL AMR structure to VTK's native AMR or multiblock data structures.
+
+1. **`vtkNonOverlappingAMR`**: Since BATL blocks do not overlap and child blocks completely replace their parents, the data maps perfectly to `vtkNonOverlappingAMR`. This format is a specialized subclass of `vtkUniformGridAMR`. Each active BATL block becomes a `vtkUniformGrid` at its respective AMR level. This is the most mathematically correct representation for BATL.
+2. **`vtkOverlappingAMR` (Berger-Colella)**: Although BATL is non-overlapping at the active leaf level, the entire tree can be represented as overlapping by including the parent blocks (which are "unused" in BATL but exist in the tree). VTK's `vtkOverlappingAMR` uses a blanking array (`IBLANK`) to hide the portions of coarse grids covered by finer grids. This might allow ParaView to use specialized dual-mesh filters that are exclusively implemented for `vtkOverlappingAMR`.
+3. **`vtkMultiBlockDataSet`**: If native AMR metadata formats (`.vth` or `.vthb`) are difficult to generate natively in Julia using current tools (e.g., `WriteVTK.jl` has limited native AMR XML support), a simpler step up from unstructured grids is to export each BATL block as an individual structured grid (`vtkImageData` or `vtkRectilinearGrid`) and group them into a `vtkMultiBlockDataSet` (`.vtm` file). This completely avoids the unstructured connectivity overhead, although it still lacks explicit AMR level metadata.
+
+To implement native VTK AMR, one would need to write a `.vth` / `.vthb` (Hierarchical Box) XML file, providing the exact bounding boxes and refinement levels corresponding to the BATL tree metadata.

src/derived.jl

@@ -274,6 +274,7 @@ end
 Return the electric field ``\\mathbf{E}_{p_e} = -\\frac{1}{n_e e} \\nabla \\cdot \\mathbf{P}_e``
 derived from the divergence of the electron pressure tensor. Units are in μV/m if
 PLANETARY or km units are used, otherwise 1.0.
+Note that species 0 is by convention electrons, but it is not guaranteed.
 """
 @inline function get_pe_E(bd::BatsrusIDL, species = 0; mass = nothing)
     return _get_pe_E(bd, species, mass)

src/select.jl

@@ -252,7 +252,8 @@ end
 end
 
 @inline function get_vectors_indices(bd::BatsrusIDL, ::Val{:J})
-    return (findindex(bd, "jx"), findindex(bd, "jy"), findindex(bd, "jz"))
+    idx = findindex(bd, "jx")
+    return (idx, idx + 1, idx + 2)
 end
 
 function get_vectors_indices(bd::BatsrusIDL, ::Val{V}) where {V}

src/vtk.jl

@@ -245,7 +245,7 @@ function readtree(filetag)
 end
 
 """
-     find_grid_block(batl, xyz_D)
+    find_grid_block(batl, xyz_D)
 
 Return processor local block index that contains a point. Input location should be given in
 Cartesian coordinates.
@@ -279,7 +279,7 @@ function find_grid_block(batl::Batl, xyz_D)
 end
 
 """
-     find_tree_node(batl, Coord_D)
+    find_tree_node(batl, Coord_D)
 
 Find the node that contains a point. The point coordinates should be given in generalized
 coordinates normalized by the domain size.
@@ -347,7 +347,7 @@ function ibits(i, pos, len)
 end
 
 """
-     order_tree(batl)
+    order_tree(batl)
 
 Return maximum AMR level in the used block and the Morton curve order.
 Set iNodeMorton_I indirect index arrays according to
@@ -391,7 +391,7 @@ function order_tree(batl::Batl)
 end
 
 """
-     order_children!(batl::Batl, iNode, iMorton::Int, iNodeMorton_I::Vector{Int32})
+    order_children!(batl::Batl, iNode, iMorton::Int, iNodeMorton_I::Vector{Int32})
 
 Recursively apply Morton ordering for nodes below a root block.
 Store result into iNodeMorton_I and iMortonNode_A using the iMorton index.