diff --git a/Project.toml b/Project.toml
index 0c531114a..adaf963bc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,6 +24,9 @@ Fortran90Namelists = "8fb689aa-71ff-4044-8071-0cffc910b57d"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FusionMaterials = "4c86da02-02c8-4634-8460-96566129f8e0"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 IMAS = "13ead8c1-b7d1-41bb-a6d0-5b8b65ed587a"
 IMASDD = "06b86afa-9f21-11ec-2ef8-e51b8960cfc5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -41,6 +44,7 @@ Metaheuristics = "bcdb8e00-2c21-11e9-3065-2b553b22f898"
 MillerExtendedHarmonic = "c82744c2-dc08-461a-8c37-87ab04d0f9b8"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 NNeutronics = "a9424c20-d414-11ec-167b-9106c24d956c"
+NeutronTransport = "161f4978-d3bb-4696-90f4-058341551711"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
@@ -52,13 +56,18 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 QED = "8bcbec86-48c7-11ec-16a6-c1bc83299373"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RayTracing = "3007c720-8091-40db-9339-e09e4eb4c7ea"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SimulationParameters = "32ab26d3-25d6-405d-a295-367385e16093"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TAUENN = "4f5ffc88-c87f-11eb-142e-4d2ba0394a67"
 TEQUILA = "a60c9cbd-e72f-4185-96b6-b8fc312c4d37"
 TGLFNN = "558c7b13-fd9f-4806-b461-592296cfa9d0"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 VacuumFields = "9d9223b5-c5da-4bf4-abee-2a8bb6775a49"
 Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9"
 
diff --git a/playground/2d_neutron_transport_example.ipynb b/playground/2d_neutron_transport_example.ipynb
new file mode 100644
index 000000000..ee8632f14
--- /dev/null
+++ b/playground/2d_neutron_transport_example.ipynb
@@ -0,0 +1,90 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1c50669",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using Revise\n",
+    "using FUSE"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "009f8417",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ini, act = FUSE.case_parameters(:FPP; version=:v1_demount, init_from=:scalars);\n",
+    "dd = FUSE.init(ini, act; do_plot=false);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "02231999",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_path=\"/fusion/projects/ird/ptp/mclaughlink/\"\n",
+    "data_filename=\"combined_mgxs_lib.h5\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f81ec0c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@time sol = FUSE.neutron_transport_2d(dd, data_path, data_filename)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a6195600",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "FUSE.get_tbr(sol, dd, data_path, data_filename)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e6d8682a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "FUSE.write_neutron_flux_vtk(sol)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "931b8a69",
+   "metadata": {},
+   "source": [
+    "You can now use Paraview or a similar program to view the fluxes.vtu file. If running on Omega, you can type `/fusion/projects/codes/neutronics/ParaView-5.11.1-MPI-Linux-Python3.9-x86_64/bin/paraview` in your terminal window to launch Paraview."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.9.2",
+   "language": "julia",
+   "name": "julia-1.9"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/playground/neutron_transport_example.ipynb b/playground/neutron_transport_example.ipynb
new file mode 100644
index 000000000..de3173126
--- /dev/null
+++ b/playground/neutron_transport_example.ipynb
@@ -0,0 +1,101 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using Revise\n",
+    "using FUSE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Build an arbitrary 1D geometry using `FUSE.concentric_cylinders()`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layer_thicknesses =  [100, 0.5,20,1]  # eventually replace with dd.build.layer.ΔR ?\n",
+    "material_names = [\"DT_plasma\", \"Tungsten\", \"lithium-lead-Li6enrich=90\",\"eurofer\"] # eventually replace with dd.build.layer.material\n",
+    "layers = [\"plasma\", \"First Wall\", \"Breeder\", \"Back wall\"]\n",
+    "mat_tags = Dict(zip(layers, material_names))\n",
+    "data_path = \"/fusion/projects/ird/ptp/mclaughlink/\"\n",
+    "data_filename= \"combined_mgxs_lib.h5\"\n",
+    "save_path = data_path\n",
+    "save_filename = \"concentric_circles\"\n",
+    "sol = FUSE.concentric_circles(layer_thicknesses, material_names, layers, data_path, data_filename, save_path, save_filename)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tbr = FUSE.get_tbr(sol, mat_tags, data_path, data_filename)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Or, to build a 1D simulation based on FUSE's midplane layer thicknesses, use `FUSE.neutron_transport_1d()`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ini, act = FUSE.case_parameters(:FPP; version=:v1_demount, init_from=:scalars);\n",
+    "dd = FUSE.init(ini, act; do_plot=false);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol = FUSE.neutron_transport_1d(dd, data_path, data_filename)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tbr = FUSE.get_tbr(sol, dd, data_path, data_filename)"
+   ]
+  }
+ ],
+ "metadata": {
+  "@webio": {
+   "lastCommId": null,
+   "lastKernelId": null
+  },
+  "kernelspec": {
+   "display_name": "Julia (10 threads) 1.9.0",
+   "language": "julia",
+   "name": "julia-_10-threads_-1.9"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/src/FUSE.jl b/src/FUSE.jl
index ccfdd030a..f8193a515 100644
--- a/src/FUSE.jl
+++ b/src/FUSE.jl
@@ -62,6 +62,7 @@ include(joinpath("actors", "build", "cx_actor.jl"))
 
 include(joinpath("actors", "nuclear", "blanket_actor.jl"))
 include(joinpath("actors", "nuclear", "neutronics_actor.jl"))
+include(joinpath("actors", "nuclear", "neutron_transport_actor.jl"))
 
 include(joinpath("actors", "current", "qed_actor.jl"))
 include(joinpath("actors", "current", "steadycurrent_actor.jl"))
diff --git a/src/actors/nuclear/neutron_transport_actor.jl b/src/actors/nuclear/neutron_transport_actor.jl
new file mode 100644
index 000000000..246ca9119
--- /dev/null
+++ b/src/actors/nuclear/neutron_transport_actor.jl
@@ -0,0 +1,420 @@
+using HDF5
+using NeutronTransport
+using GridapGmsh
+using GridapGmsh: gmsh, GmshDiscreteModel
+using Gridap
+import NeutronTransport: MoCSolution
+import Gridap: writevtk
+import Gridap.Geometry: get_triangulation
+
+const XSs = CrossSections
+
+function get_xss_from_hdf5(
+    filename::String,
+    material::String,
+    material_tag::String=material;
+    )
+    file = h5open(filename, "r")
+    νΣf = read(file["material"][material]["nu-fission"]["average"])
+    Σt = read(file["material"][material]["nu-transport"]["average"])
+    Σs0 = permutedims(read(file["material"][material]["consistent nu-scatter matrix"]["average"]))
+    return CrossSections(material_tag, length(νΣf); νΣf, Σt, Σs0)
+end
+
+function get_xss_from_hdf5(filename::String, material::String, xs_type::Vector{String})
+    xs_dict = Dict{String,Vector{Float64}}()
+    file = h5open(filename, "r")
+    for xs_name in xs_type
+        xs = read(file["material"][material][xs_name]["average"])
+        xs_dict[xs_name] = xs
+    end
+    close(file)
+    return xs_dict
+end
+
+function get_xss_from_hdf5(filename::String, materials::Vector{String}, xs_type::Vector{String})
+    xs_dict=Dict{String, Dict{String,Vector{Float64}}}()
+    for mat in materials
+        mat_xs_dict = get_xss_from_hdf5(filename, mat, xs_type)
+        xs_dict[mat]=mat_xs_dict
+    end
+    return xs_dict
+end
+
+function concentric_circles(layer_thicknesses::Vector{T}, materials::Vector{String}, layers::Vector{String}, data_path::String, data_filename::String, save_path::String=pwd(), save_filename::String="concentric_circles") where T
+    gmsh.initialize()
+    model_bound = sum(layer_thicknesses) + 1
+
+    mats = deepcopy(materials)
+    thicknesses = deepcopy(layer_thicknesses)
+    names = deepcopy(layers)
+    push!(mats, "Air-(dry,-near-sea-level)")
+    push!(names, "Exterior")
+
+    append!(thicknesses, 1.0)
+    factory = gmsh.model.geo
+    lc = model_bound/30
+
+    radius=0.0
+
+    data_filename=joinpath(data_path, data_filename)
+
+    for (idx, thickness) in enumerate(thicknesses)
+        radius+=thickness
+        mats[idx] = replace(mats[idx], " " => "-")
+        if mats[idx] == "Vacuum"
+            mats[idx] = "Air-(dry,-near-sea-level)"
+        end 
+        material = factory.addPhysicalGroup(2, [idx], idx)
+        GridapGmsh.gmsh.model.setPhysicalName(2, idx, names[idx])
+
+        if idx == length(mats)
+            # square cell
+            top_right = factory.addPoint(radius, radius, 0, lc, 5*(idx-1)+1)
+            bottom_right = factory.addPoint(radius, -radius, 0, lc, 5*(idx-1)+2)
+            bottom_left = factory.addPoint(-radius, -radius, 0, lc, 5*(idx-1)+3)
+            top_left = factory.addPoint(-radius, radius, 0, lc, 5*(idx-1)+4)
+            right = factory.addLine(top_right, bottom_right,  4*(idx-1)+1)
+            bottom = factory.addLine(bottom_right, bottom_left, 4*(idx-1)+2)
+            left = factory.addLine(bottom_left, top_left, 4*(idx-1)+3)
+            top = factory.addLine(top_left, top_right, 4*(idx-1)+4)
+            boundary = factory.addCurveLoop([right, bottom, left, top], idx)
+
+            # boundaries
+            right_bound = factory.addPhysicalGroup(1,  [right], material+1)
+            bottom_bound = factory.addPhysicalGroup(1, [bottom], material+2)
+            left_bound = factory.addPhysicalGroup(1, [left], material+3)
+            top_bound = factory.addPhysicalGroup(1, [top], material+4)
+            GridapGmsh.gmsh.model.setPhysicalName(1, right_bound, "right")
+            GridapGmsh.gmsh.model.setPhysicalName(1, bottom_bound, "bottom")
+            GridapGmsh.gmsh.model.setPhysicalName(1, left_bound, "left")
+            GridapGmsh.gmsh.model.setPhysicalName(1, top_bound, "top")
+        else
+            # make inner circle points
+            center = factory.addPoint(0, 0, 0, lc, 5*(idx-1)+1)
+            right = factory.addPoint(radius, 0, 0, lc, 5*(idx-1)+2)
+            top = factory.addPoint(0, radius, 0, lc, 5*(idx-1)+3)
+            left = factory.addPoint(-radius, 0, 0, lc, 5*(idx-1)+4)
+            bottom = factory.addPoint(0, -radius, 0, lc, 5*(idx-1)+5)
+
+            # make arcs and circle
+            right_top = factory.addCircleArc(right, center, top, 4*(idx-1)+1)
+            top_left = factory.addCircleArc(top, center, left, 4*(idx-1)+2)
+            left_bottom = factory.addCircleArc(left, center, bottom, 4*(idx-1)+3)
+            bottom_right = factory.addCircleArc(bottom, center, right, 4*(idx-1)+4)
+            circle = factory.addCurveLoop([right_top, top_left, left_bottom, bottom_right], idx)
+        end
+
+        # make surfaces and materials
+        if idx==1
+            surface = factory.addPlaneSurface([idx], idx)
+        else
+            surface = factory.addPlaneSurface([idx, idx-1], idx)
+        end
+    end
+
+    xss = [get_xss_from_hdf5(data_filename, mats[idx], names[idx]) for idx in eachindex(thicknesses)]
+
+    factory.synchronize()
+
+    gmsh.model.mesh.generate(2)
+
+    mshfile = joinpath(save_path, save_filename * ".msh")
+
+    gmsh.write(mshfile)
+
+    # gmsh.fltk.run()
+
+    gmsh.finalize()
+
+    model = GmshDiscreteModel(mshfile; renumber=true)
+    jsonfile = joinpath(save_path, save_filename * ".json")
+    Gridap.Io.to_json_file(model, jsonfile)
+
+    jsonfile = joinpath(save_path, save_filename * ".json")
+    geometry = DiscreteModelFromFile(jsonfile)
+
+    # number of azimuthal angles
+    nφ = 4
+
+    # azimuthal spacing
+    δ = 3.0 
+
+    # boundary conditions
+    bcs = BoundaryConditions(top=Vaccum, left=Vaccum, bottom=Vaccum, right=Vaccum)
+
+    # initialize track generator
+    tg = TrackGenerator(geometry, nφ, δ, bcs=bcs)
+    
+    # perform ray tracing
+    trace!(tg)
+
+    # proceed to segmentation
+    segmentize!(tg)
+
+    # polar quadrature
+    pq = NeutronTransport.TabuchiYamamoto(2)
+
+    # define the problem
+    prob = MoCProblem(tg, pq, xss)
+
+    # define fixed source material
+    fixed_sources = set_fixed_source_material(prob, "plasma", 8 , 1)
+
+    # solve
+    sol = NeutronTransport.solve(prob, fixed_sources, debug=true, max_residual=0.05, max_iterations=500)
+
+    return sol
+end
+
+"""
+    get_cell_φ(sol, cell, g)
+
+Returns the average φ in a given cell and energy group.
+"""
+function get_cell_φ(sol::MoCSolution{T}, cell::Int, g::Int) where T
+    φs = sol(g)[findall(c -> c == cell, sol.prob.fsr_tag)]
+    # vols = sol.prob.trackgenerator.volumes[findall(c -> c == cell, sol.prob.fsr_tag)]
+    
+    # vol_integrated_φ = sum(φs .* vols / 31415)
+    # return vol_integrated_φ
+
+    mean_φ = sum(φs)/length(φs)
+    return mean_φ
+end
+
+function get_cell_φ(sol::MoCSolution{T}, cell::Int, gs::UnitRange{Int64}) where T
+    φs = Float64[]
+    for g in gs
+        append!(φs, get_cell_φ(sol, cell, g))
+    end
+    return φs
+end
+
+function get_tbr(sol::MoCSolution{T}, mat_tags::Dict{String, String}, data_path::String, data_filename::String) where T
+    NeutronTransport.@unpack φ, prob = sol
+    NeutronTransport.@unpack trackgenerator, fsr_tag, xss = prob
+    NeutronTransport.@unpack volumes = trackgenerator
+    NGroups = NeutronTransport.ngroups(prob)
+    NRegions = NeutronTransport.nregions(prob)
+
+    vol_integrated_φ = Vector{Float64}()
+    tbr=0.0
+    plasma_vol = sum(volumes[findall(t -> t == 1, fsr_tag)])
+
+    for i in 1:NRegions
+        if xss[fsr_tag[i]].name == "Exterior"
+            material = "Air-(dry,-near-sea-level)"
+        else
+            material = mat_tags[xss[fsr_tag[i]].name]
+            material = replace(material," " => "-", "Vacuum" => "Air-(dry,-near-sea-level)")
+        end
+        tbr_xs_dict = get_xss_from_hdf5(joinpath(data_path,data_filename), material, ["(n,Xt)"])
+        for g in 1:NGroups
+            ig = NeutronTransport.@region_index(i, g)
+            push!(vol_integrated_φ, sol.φ[ig] * volumes[i] / plasma_vol)
+            tbr+=vol_integrated_φ[ig]*tbr_xs_dict["(n,Xt)"][g]
+        end        
+    end
+    return tbr
+end
+
+function get_tbr(sol::MoCSolution{T}, dd::IMAS.dd, data_path::String, data_filename::String) where T
+    NeutronTransport.@unpack φ, prob = sol
+    NeutronTransport.@unpack trackgenerator, fsr_tag, xss = prob
+    NeutronTransport.@unpack volumes = trackgenerator
+    NGroups = NeutronTransport.ngroups(prob)
+    NRegions = NeutronTransport.nregions(prob)
+    
+    vol_integrated_φ = Vector{Float64}()
+    tbr=0.0
+    plasma_vol = sum(volumes[findall(t -> t == 1, fsr_tag)])
+
+    for i in 1:NRegions
+        if xss[fsr_tag[i]].name == "Exterior"
+            material = "Air-(dry,-near-sea-level)"
+        else
+            material = IMAS.get_build_layers(dd.build.layer, name=xss[fsr_tag[i]].name)[1].material
+            material = replace(material," " => "-", "Vacuum" => "Air-(dry,-near-sea-level)")
+        end
+        tbr_xs_dict = get_xss_from_hdf5(joinpath(data_path,data_filename), material, ["(n,Xt)"])
+        for g in 1:NGroups
+            ig = NeutronTransport.@region_index(i, g)
+            push!(vol_integrated_φ, sol.φ[ig] * volumes[i] / plasma_vol)
+            tbr+=vol_integrated_φ[ig]*tbr_xs_dict["(n,Xt)"][g]
+        end        
+    end
+    return tbr
+end
+    
+
+function write_neutron_flux_vtk(sol::MoCSolution{T}, groups::Union{Vector{Int},UnitRange{Int}}=1:length(sol.prob.xss[1].Σt)) where T
+    NeutronTransport.@unpack prob = sol
+    NeutronTransport.@unpack trackgenerator = prob
+    triang = get_triangulation(trackgenerator.mesh.model)
+    writevtk(triang, "fluxes", cellfields=[string(g) => sol(g) for g in groups])
+end
+
+function neutron_transport_1d(dd::IMAS.dd, data_path::String, data_filename::String, save_path::String=pwd(), save_filename::String="fuse_neutron_transport")
+    thicknesses = [layer.thickness*100 for layer in IMAS.get_build_layers(dd.build.layer, fs=[IMAS._in_, IMAS._hfs_, IMAS._lhfs_])]
+    materials = [layer.material for layer in IMAS.get_build_layers(dd.build.layer, fs=[IMAS._in_, IMAS._hfs_, IMAS._lhfs_])]
+    layers = [layer.name for layer in IMAS.get_build_layers(dd.build.layer, fs=[IMAS._in_, IMAS._hfs_, IMAS._lhfs_])]
+    return concentric_circles(thicknesses, materials, layers, data_path, data_filename, save_path, save_filename)
+end
+
+function construct_boundary(r_coords, z_coords, surface_number)
+    factory = gmsh.model.geo
+    lc = 30
+    coords = [factory.addPoint(r,z,0,lc,surface_number*100000+x) for (x, (r, z)) in enumerate(zip(r_coords, z_coords))]
+    pop!(coords)
+    lines = [factory.addLine(coords[x], coords[x+1], surface_number*100000+x) for x in 1:length(coords)-1]
+    append!(lines, factory.addLine(coords[length(coords)], coords[1], surface_number*100000+length(coords)))
+    boundary = factory.addCurveLoop(lines, surface_number)
+    return boundary
+end
+
+function construct_boundary(r_coords, z_coords)
+    # square cell
+    factory = gmsh.model.geo
+    lc = 30
+    top_right = factory.addPoint(r_coords[1], z_coords[1], 0, lc, 999991)
+    bottom_right = factory.addPoint(r_coords[2], z_coords[2], 0, lc, 999992)
+    bottom_left = factory.addPoint(r_coords[3], z_coords[3], 0, lc, 999993)
+    top_left = factory.addPoint(r_coords[4], z_coords[4], 0, lc, 999994)
+    right = factory.addLine(top_right, bottom_right,  999995)
+    bottom = factory.addLine(bottom_right, bottom_left, 999996)
+    left = factory.addLine(bottom_left, top_left, 999997)
+    top = factory.addLine(top_left, top_right, 999998)
+    boundary = factory.addCurveLoop([right, bottom, left, top], 999999)
+
+    # boundaries
+    right_bound = factory.addPhysicalGroup(1,  [right], 999991)
+    bottom_bound = factory.addPhysicalGroup(1, [bottom], 999992)
+    left_bound = factory.addPhysicalGroup(1, [left], 999993)
+    top_bound = factory.addPhysicalGroup(1, [top], 999994)
+    GridapGmsh.gmsh.model.setPhysicalName(1, right_bound, "right")
+    GridapGmsh.gmsh.model.setPhysicalName(1, bottom_bound, "bottom")
+    GridapGmsh.gmsh.model.setPhysicalName(1, left_bound, "left")
+    GridapGmsh.gmsh.model.setPhysicalName(1, top_bound, "top")
+    return boundary
+end
+
+function construct_circle(center_coords::Vector, radius, surface_id)
+    # make inner circle points
+    factory = gmsh.model.geo
+    lc = 30
+    center = factory.addPoint(center_coords[1], center_coords[2], 0, lc, 5*(surface_id-1)+1)
+    right = factory.addPoint(center_coords[1]+radius, center_coords[2], 0, lc, 5*(surface_id-1)+2)
+    top = factory.addPoint(center_coords[1], center_coords[2]+radius, 0, lc, 5*(surface_id-1)+3)
+    left = factory.addPoint(center_coords[1]-radius, center_coords[2], 0, lc, 5*(surface_id-1)+4)
+    bottom = factory.addPoint(center_coords[1], center_coords[2]-radius, 0, lc, 5*(surface_id-1)+5)
+
+    # make arcs and circle
+    right_top = factory.addCircleArc(right, center, top, 4*(surface_id-1)+1)
+    top_left = factory.addCircleArc(top, center, left, 4*(surface_id-1)+2)
+    left_bottom = factory.addCircleArc(left, center, bottom, 4*(surface_id-1)+3)
+    bottom_right = factory.addCircleArc(bottom, center, right, 4*(surface_id-1)+4)
+    circle = factory.addCurveLoop([right_top, top_left, left_bottom, bottom_right], surface_id)
+    return circle
+end
+
+function neutron_transport_2d(dd::IMAS.dd, data_path::String, data_filename::String, save_path::String=pwd(), save_filename::String="fuse_neutron_transport")
+    gmsh.initialize()
+    factory = gmsh.model.geo
+    layers = IMAS.get_build_layers(dd.build.layer, fs=[IMAS._hfs_, IMAS._lhfs_]) 
+    pushfirst!(layers, IMAS.get_build_layer(dd.build.layer, type=IMAS._oh_))
+    layers = [layer for layer in layers if !occursin("first", layer.name)] # temporary fix
+    names = [replace(layer.name, "plasma" => "plsm") for layer in layers]
+    materials = [replace(layer.material," " => "-", "Vacuum" => "Air-(dry,-near-sea-level)") for layer in layers]
+    boundaries = [construct_boundary(layer.outline.r * 100, layer.outline.z * 100, idx) for (idx, layer) in enumerate(layers)]
+    equilibrium_magnetic_axis = dd.equilibrium.time_slice[2].global_quantities.magnetic_axis
+    push!(names, "plasma")
+    push!(materials, "DT_plasma")
+    push!(boundaries, construct_circle([equilibrium_magnetic_axis.r * 100, equilibrium_magnetic_axis.z * 100], 50, 88888))
+    push!(materials, "Air-(dry,-near-sea-level)")
+    push!(names, "Exterior")
+    max_z = maximum([maximum(layer.outline.z) for layer in layers]) * 100 + 1.0
+    min_z = minimum([minimum(layer.outline.z) for layer in layers]) * 100 - 1.0
+    max_r = maximum([maximum(layer.outline.r) for layer in layers]) * 100 + 1.0
+    min_r = 0.0
+    r_bounds = [max_r, min_r, min_r, max_r]
+    z_bounds = [max_z, max_z, min_z, min_z]
+    push!(boundaries, construct_boundary(r_bounds, z_bounds))
+    outer_bounds = [last(boundaries)]
+    counter = 1
+
+    for idx in eachindex(boundaries)
+        if idx <= length(layers) && (layers[idx].fs == Int(IMAS._hfs_) || layers[idx].fs == Int(IMAS._lhfs_))
+            surface = factory.addPlaneSurface([boundaries[idx], boundaries[idx+1]], idx)
+            if counter == 1
+                push!(outer_bounds, surface)
+            end
+        elseif idx == length(boundaries)
+            surface = factory.addPlaneSurface(outer_bounds, idx)
+        else
+            surface = factory.addPlaneSurface([boundaries[idx]], idx)
+            if idx <= length(layers) && layers[idx].fs == Int(IMAS._in_)
+                push!(outer_bounds, surface)
+            end
+        end
+        material = factory.addPhysicalGroup(2, [idx], idx)
+        GridapGmsh.gmsh.model.setPhysicalName(2, idx, names[idx]) 
+    end
+
+    data_filename=joinpath(data_path, data_filename)
+    xss = [get_xss_from_hdf5(data_filename, materials[idx], names[idx]) for idx in eachindex(boundaries)]
+    factory.synchronize()
+    gmsh.model.mesh.generate(2)
+    mshfile = joinpath(save_path, save_filename * ".msh")
+    gmsh.fltk.run()
+    gmsh.write(mshfile)
+    gmsh.finalize()
+    
+    sol = run_2d(dd, xss, mshfile, save_path, save_filename)
+    return sol
+end
+
+function run_2d(dd, xss, mshfile, save_path::String=pwd(), save_filename::String="fuse_neutron_transport")
+    model = GmshDiscreteModel(mshfile; renumber=true)
+    jsonfile = joinpath(save_path, save_filename * ".json")
+    Gridap.Io.to_json_file(model, jsonfile)
+
+    jsonfile = joinpath(save_path, save_filename * ".json")
+    geometry = DiscreteModelFromFile(jsonfile)
+
+    # number of azimuthal angles
+    nφ = 16
+
+    # azimuthal spacing
+    δ = 1.5 
+
+    # boundary conditions
+    bcs = BoundaryConditions(top=Vaccum, left=Vaccum, bottom=Vaccum, right=Vaccum)
+
+    # initialize track generator
+    tg = TrackGenerator(geometry, nφ, δ, bcs=bcs)
+    
+    # perform ray tracing
+    trace!(tg)
+
+    # proceed to segmentation
+    segmentize!(tg)
+
+    # polar quadrature
+    pq = NeutronTransport.TabuchiYamamoto(2)
+
+    # define the problem
+    prob = MoCProblem(tg, pq, xss)
+
+    # define fixed source material
+    if IMAS.get_build_layers(dd.build.layer, type=_plasma_)[1].material == "DD_plasma"
+        fixed_sources = set_fixed_source_material(prob, "plasma", 42 , 1)
+    else
+        fixed_sources = set_fixed_source_material(prob, "plasma", 8 , 1)
+    end
+
+    # solve
+    sol = NeutronTransport.solve(prob, fixed_sources, debug=true, max_residual=0.05, max_iterations=500)
+
+    return sol
+end
\ No newline at end of file