15 improve atlas

README.md

@@ -14,10 +14,10 @@ export GRIDTOOLS_JL_PATH="..."
 export GT4PY_PATH="..."
 # create python virtual environemnt
 #  make sure to use a python version that is compatible with GT4Py
-python -m venv venv
+python -m venv .venv
 # activate virtual env
 #  this command has be run everytime GridTools.jl is used
-source venv/bin/activate
+source .venv/bin/activate
 # clone gt4py
 git clone --branch fix_python_interp_path_in_cmake [email protected]:tehrengruber/gt4py.git
 #git clone [email protected]:GridTools/gt4py.git $GT4PY_PATH

advection/README.md

@@ -0,0 +1,48 @@
+### README for Running `advection_example.jl`
+
+This README provides instructions on how to run the `advection_example.jl` script for simulating advection using the Atlas library. The script allows for terminal visualization, which can be enabled as described below.
+
+#### Prerequisites
+
+1. **Python Environment and Atlas4py Installation**:
+   - Ensure that your Python environment is activated.
+     ```sh
+     source .venv/bin/activate # Ignore if the env is already activated
+     ```
+   - Install the `atlas4py` package using the following command:
+     ```sh
+     pip install -i https://test.pypi.org/simple/ atlas4py
+     ```
+
+2. **Enabling Visualization** (optional):
+   - The script has a `VISUALIZATION_FLAG` that can be set to enable or disable visualization on the terminal. Ensure that this flag is set to `true` in the `advection_example.jl` script if you wish to enable visualization.
+   - Note: Other parameters such as the number of iterations can be changed in the `# Simulation Parameters` section of the script.
+
+#### Running the Simulation
+
+1. **Running the Script**:
+   - Use the following command to run the `advection_example.jl` script with Julia:
+     ```sh
+     julia --color=yes --project=$GRIDTOOLS_JL_PATH/GridTools.jl $GRIDTOOLS_JL_PATH/GridTools.jl/src/examples/advection/advection_example.jl
+     ```
+
+#### Example
+
+Here is an example of how to set the `VISUALIZATION_FLAG` in the `advection_example.jl` script and run the simulation:
+
+1. **Setting the Visualization Flag**:
+   - Open the `advection_example.jl` script.
+   - Set the `VISUALIZATION_FLAG` to `true`:
+     ```julia
+     const VISUALIZATION_FLAG = true
+     ```
+
+2. **Running the Simulation**:
+   - Open your terminal.
+   - Run the script with the following command:
+     ```sh
+     export GRIDTOOLS_JL_PATH=...
+     julia --color=yes --project=. $GRIDTOOLS_JL_PATH/src/examples/advection/advection_example.jl
+     ```
+
+By following these steps, you should be able to run the `advection_example.jl` script and visualize the advection simulation results on your terminal.

src/atlas/advection_test.jl → advection/advection_miniapp.jl

@@ -1,17 +1,25 @@
-# advection test
+# Advection Example
+# This script demonstrates an advection simulation using the Atlas library.
+
 using Printf
 using Debugger
 using Statistics
 using Profile
 using GridTools
 
+const global VISUALIZATION_FLAG::Bool=false
+
+# Mesh Definitions --------------------------------------------------------------------------------------------
+# Define dimensions for the mesh
 Cell_ = Dimension{:Cell_, HORIZONTAL}
 Edge_ = Dimension{:Edge_, HORIZONTAL}
 Vertex_ = Dimension{:Vertex_, HORIZONTAL}
 K_ = Dimension{:K_, VERTICAL}
 V2VDim_ = Dimension{:V2V_, LOCAL}
 V2EDim_ = Dimension{:V2E_, LOCAL}
 E2VDim_ = Dimension{:E2V_, LOCAL}
+
+# Instantiate dimension objects
 Cell = Cell_()
 K = K_()
 Edge = Edge_()
@@ -20,41 +28,53 @@ V2VDim = V2VDim_()
 V2EDim = V2EDim_()
 E2VDim = E2VDim_()
 
+# Define field offsets to describe the relationships between different dimensions
 V2V = FieldOffset("V2V", source = Vertex, target = (Vertex, V2VDim))
 E2V = FieldOffset("E2V", source = Vertex, target = (Edge, E2VDim))
 V2E = FieldOffset("V2E", source = Edge, target = (Vertex, V2EDim))
 Koff = FieldOffset("Koff", source = K, target = K)
 
-include("atlas_mesh.jl")
+# Include additional necessary files for mesh, state container, metric calculations, and advection operations
+include("../src/atlas/atlas_mesh.jl")
 include("state_container.jl")
 include("metric.jl")
 include("advection.jl")
+include("visualization_utils.jl")
 
-
+# Grid and Mesh Initialization --------------------------------------------------------------------------------
+# Create a structured grid and mesh for the simulation
 grid = atlas.StructuredGrid("O50")
 mesh = AtlasMesh(grid, num_level = 30)
 
+# Simulation Parameters ---------------------------------------------------------------------------------------
 δt = 1800.0  # time step in s
 niter = 50
 eps = 1.0e-8
 
+# Calculate metric properties from the mesh
 metric = m_from_mesh(mesh)
+
+# Define the spatial extent of the mesh
 origin = minimum(mesh.xyarc, dims = 1)
 extent = maximum(mesh.xyarc, dims = 1) .- minimum(mesh.xyarc, dims = 1)
 xlim = (minimum(mesh.xyarc[:, 1]), maximum(mesh.xyarc[:, 1]))
 ylim = (minimum(mesh.xyarc[:, 2]), maximum(mesh.xyarc[:, 2]))
 
+# Get dimensions of various elements in the mesh
 vertex_dim = getproperty(mesh, DIMENSION_TO_SIZE_ATTR[Vertex])
 k_dim = getproperty(mesh, DIMENSION_TO_SIZE_ATTR[K])
 edge_dim = getproperty(mesh, DIMENSION_TO_SIZE_ATTR[Edge])
 
+# Define vertical level indices
 level_indices = Field(K, zeros(Int, k_dim))
 level_indices .= collect(0:mesh.num_level-1)
 
-
+# Initialize state containers for the current and next states
 state = sc_from_mesh(mesh)
 state_next = sc_from_mesh(mesh)
 
+# Temporary Fields Initialization -----------------------------------------------------------------------------
+# Create temporary fields used in the computation
 tmp_fields = Dict{String, Field}()
 for i = 1:6
     tmp_fields[@sprintf("tmp_vertex_%d", i)] = Field((Vertex, K), zeros(vertex_dim, k_dim))
@@ -63,33 +83,40 @@ for j = 1:3
     tmp_fields[@sprintf("tmp_edge_%d", j)] = Field((Edge, K), zeros(edge_dim, k_dim))
 end
 
+# Initial Conditions -------------------------------------------------------------------------------------------
+# Define the initial conditions for the scalar field (rho) using a field operator
 @field_operator function initial_rho(
     mesh_radius::Float64,
     mesh_xydeg_x::Field{Tuple{Vertex_}, Float64},
     mesh_xydeg_y::Field{Tuple{Vertex_}, Float64},
     mesh_vertex_ghost_mask::Field{Tuple{Vertex_}, Bool}
 )::Field{Tuple{Vertex_, K_}, Float64}
+    # Define constants for the initial condition
     lonc = 0.5 * pi
     latc = 0.0
     _deg2rad = 2.0 * pi / 360.0
 
+    # Convert mesh coordinates from degrees to radians
     mesh_xyrad_x, mesh_xyrad_y = mesh_xydeg_x .* _deg2rad, mesh_xydeg_y .* _deg2rad
     rsina, rcosa = sin.(mesh_xyrad_y), cos.(mesh_xyrad_y)
 
+    # Compute the distance from the center point (lonc, latc)
     zdist =
         mesh_radius .*
         acos.(sin(latc) .* rsina .+ cos(latc) .* rcosa .* cos.(mesh_xyrad_x .- lonc))
 
+    # Calculate the radial profile
     rpr = (zdist ./ (mesh_radius / 2.0)) .^ 2.0
     rpr = min.(1.0, rpr)
 
+    # Return the initial scalar field values, setting ghost cells to zero
     return broadcast(
         where(mesh_vertex_ghost_mask, 0.0, 0.5 .* (1.0 .+ cos.(pi .* rpr))),
         (Vertex, K)
     )
 end
 
-
+# Initialize the scalar field (rho) using the initial_rho function
 initial_rho(
     mesh.radius,
     mesh.xydeg_x,
@@ -99,6 +126,7 @@ initial_rho(
     offset_provider = mesh.offset_provider
 )
 
+# Define the initial conditions for the velocity field using a field operator
 @field_operator function initial_velocity(
     mesh_xydeg_x::Field{Tuple{Vertex_}, Float64},
     mesh_xydeg_y::Field{Tuple{Vertex_}, Float64},
@@ -110,22 +138,31 @@ initial_rho(
     Field{Tuple{Vertex_, K_}, Float64},
     Field{Tuple{Vertex_, K_}, Float64}
 }
+    # Convert mesh coordinates from degrees to radians
     _deg2rad = 2.0 * pi / 360.0
     mesh_xyrad_x, mesh_xyrad_y = mesh_xydeg_x .* _deg2rad, mesh_xydeg_y .* _deg2rad
+
+    # Set initial velocity parameters
     u0 = 22.238985328911745
     flow_angle = 0.0 * _deg2rad  # radians
 
+    # Calculate sine and cosine of mesh coordinates and flow angle
     rsina, rcosa = sin.(mesh_xyrad_y), cos.(mesh_xyrad_y)
     cosb, sinb = cos(flow_angle), sin(flow_angle)
+
+    # Compute velocity components
     uvel_x = u0 .* (cosb .* rcosa .+ rsina .* cos.(mesh_xyrad_x) .* sinb)
     uvel_y = -u0 .* sin.(mesh_xyrad_x) .* sinb
 
+    # Broadcast velocity components across the mesh
     vel_x = broadcast(uvel_x .* metric_g11 .* metric_gac, (Vertex, K))
     vel_y = broadcast(uvel_y .* metric_g22 .* metric_gac, (Vertex, K))
     vel_z = broadcast(0.0, (Vertex, K))
+
     return vel_x, vel_y, vel_z
 end
 
+# Initialize the velocity field
 initial_velocity(
     mesh.xydeg_x,
     mesh.xydeg_y,
@@ -136,10 +173,13 @@ initial_velocity(
     offset_provider = mesh.offset_provider,
 )
 
+# Copy the initial velocity field to the next state
 copyfield!(state_next.vel, state.vel)
 
+# Example of printing initial rho statistics (commented out)
 # println("min max avg of initial rho = $(minimum(state.rho.data)) , $(maximum(state.rho.data)) , $(mean(state.rho.data))")
 
+# Initialize a temporary field with vertical levels for use in the simulation
 tmp_fields["tmp_vertex_1"] .= reshape(collect(0.0:mesh.num_level-1), (1, mesh.num_level))
 nabla_z(
     tmp_fields["tmp_vertex_1"],
@@ -149,8 +189,15 @@ nabla_z(
     offset_provider = mesh.offset_provider
 )
 
-for i = 1:niter
+if VISUALIZATION_FLAG
+    # Precompute the mapping between the unstructured domain to the structured one for ASCII art visualization
+    grid_size = 50
+    mapping = precompute_mapping(mesh, xlim, ylim, grid_size)
+end
 
+# Main Simulation Loop ----------------------------------------------------------------------------------------
+for i = 1:niter
+    # Perform the upwind advection scheme to update the scalar field (rho)
     upwind_scheme(
         state.rho,
         δt,
@@ -167,16 +214,25 @@ for i = 1:niter
         offset_provider = mesh.offset_provider
     )
 
+    # Print the current timestep
     println("Timestep $i")
 
+    if VISUALIZATION_FLAG
+        # Print the current state as ASCII art every 5 timesteps
+        print_state_ascii(state, mesh, mapping, i, grid_size)
+    end
+
     # TODO: make a function out of this switch
+    # Swap the current and next state
     temp = state
     global state = state_next
     global state_next = temp
 
+    # Update the periodic boundary layers
     update_periodic_layers(mesh, state.rho)
 end
 
+# Output the final statistics for the scalar field (rho) and velocity fields
 println(
     "min max sum of final rho = $(minimum(state.rho.data)) , $(maximum(state.rho.data)) , $(sum(state.rho.data))"
 )

advection/visualization_utils.jl

@@ -0,0 +1,95 @@
+
+"""
+    precompute_mapping(mesh, xlim, ylim, grid_size)
+
+Precomputes the nearest vertex mapping for a structured grid.
+
+# Arguments
+- `mesh::Mesh`: An object containing the unstructured mesh data. It should have a property `xyarc` which is a matrix where each row represents the coordinates of a vertex.
+- `xlim::Tuple{Float64, Float64}`: A tuple `(x_min, x_max)` defining the range of x coordinates for the structured grid.
+- `ylim::Tuple{Float64, Float64}`: A tuple `(y_min, y_max)` defining the range of y coordinates for the structured grid.
+- `grid_size::Int`: The size of the structured grid (number of points along each axis).
+
+# Returns
+- `mapping::Matrix{Int}`: A matrix of size `(grid_size, grid_size)` where each element contains the index of the nearest vertex in the unstructured mesh for the corresponding point on the structured grid.
+"""
+function precompute_mapping(mesh, xlim, ylim, grid_size)
+    x_range = range(xlim[1], stop=xlim[2], length=grid_size)
+    y_range = range(ylim[1], stop=ylim[2], length=grid_size)
+    mapping = fill(0, grid_size, grid_size)
+
+    for i in 1:grid_size
+        for j in 1:grid_size
+            x = x_range[i]
+            y = y_range[j]
+            # Find the nearest vertex in the unstructured mesh
+            distances = [(mesh.xyarc[v, 1] - x)^2 + (mesh.xyarc[v, 2] - y)^2 for v in 1:size(mesh.xyarc, 1)]
+            nearest_vertex = argmin(distances)
+            mapping[i, j] = nearest_vertex
+        end
+    end
+
+    return mapping
+end
+
+"""
+    matrix_to_ascii(matrix::Matrix{Float64})::String
+
+Converts a matrix of `Float64` values to an ASCII art representation.
+
+# Arguments
+- `matrix::Matrix{Float64}`: A 2D array of `Float64` values representing the data to be converted to ASCII art.
+
+# Returns
+- `ascii_art::String`: A string containing the ASCII art representation of the matrix.
+"""
+function matrix_to_ascii(matrix::Matrix{Float64})
+    ascii_art = ""
+    chars = [' ', '.', ':', '-', '=', '+', '*', '#', '%', '@']
+    min_val = minimum(matrix)
+    max_val = maximum(matrix)
+    range_val = max_val - min_val
+
+    for row in eachrow(matrix)
+        for value in row
+            index = Int(floor((value - min_val) / range_val * (length(chars) - 1))) + 1
+            ascii_art *= chars[index]
+        end
+        ascii_art *= '\n'
+    end
+
+    return ascii_art
+end
+
+"""
+    print_state_ascii(state, mesh, mapping, timestep, grid_size=50)
+
+Prints the current state of a simulation as ASCII art.
+
+# Arguments
+- `state`: An object containing the simulation state. It should have a property `rho` with a nested property `data` which is an array of values representing the state.
+- `mesh::Mesh`: An object containing the unstructured mesh data. Used for mapping the state to grid points.
+- `mapping::Matrix{Int}`: A matrix mapping structured grid points to the nearest vertex in the unstructured mesh, as computed by the `precompute_mapping` function.
+- `timestep::Int`: An integer representing the current time step of the simulation.
+- `grid_size::Int`: (Optional) The size of the structured grid. Default is 50.
+
+# Returns
+- `Nothing`
+- This function clears the terminal and prints the ASCII art representation of the current state.
+"""
+function print_state_ascii(state, mesh, mapping, timestep, grid_size=50)
+    # Clear the terminal
+    print("\033c")
+
+    println("Timestep $timestep")
+    grid_data = zeros(Float64, grid_size, grid_size)
+
+    for i in 1:grid_size
+        for j in 1:grid_size
+            grid_data[i, j] = state.rho.data[mapping[i, j]]
+        end
+    end
+
+    ascii_art = matrix_to_ascii(grid_data)
+    println(ascii_art)
+end

docs/make.jl

@@ -17,7 +17,7 @@ makedocs(
 
 
 deploydocs(
-    repo = "github.com/jeffzwe/GridTools.jl.git",
+    repo = "github.com/GridTools/GridTools.jl.git",
     devbranch = "main",
 )