Integrate NEO with ActorNeoclassical

Makefile

@@ -30,7 +30,7 @@ CURRENTDIR := $(shell (pwd -P))
 TODAY := $(shell date +'%Y-%m-%d')
 export JULIA_NUM_THREADS ?= $(shell julia -e "println(length(Sys.cpu_info()))")
 
-FUSE_PACKAGES_MAKEFILE := ADAS BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite Fortran90Namelists FusionMaterials IMAS IMASDD MXHEquilibrium MeshTools MillerExtendedHarmonic NNeutronics QED SimulationParameters TAUENN TEQUILA TGLFNN VacuumFields 
+FUSE_PACKAGES_MAKEFILE := ADAS BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite Fortran90Namelists FusionMaterials IMAS IMASDD MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED SimulationParameters TAUENN TEQUILA TGLFNN VacuumFields 
 FUSE_PACKAGES_MAKEFILE := $(sort $(FUSE_PACKAGES_MAKEFILE))
 FUSE_PACKAGES := $(shell echo '$(FUSE_PACKAGES_MAKEFILE)' | awk '{printf("[\"%s\"", $$1); for (i=2; i<=NF; i++) printf(", \"%s\"", $$i); print "]"}')
 DEV_PACKAGES := $(shell find ../*/.git/config -exec grep ProjectTorreyPines \{\} \; | cut -d'/' -f 2 | cut -d'.' -f 1 | tr '\n' ' ')
@@ -271,6 +271,9 @@ SimulationParameters:
 BoundaryPlasmaModels:
 	$(call clone_pull_repo,$@)
 
+NEO:
+	$(call clone_pull_repo,$@)
+
 # Install IJulia
 IJulia:
 	julia -e '\

Project.toml

@@ -39,6 +39,7 @@ Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
 MeshTools = "804828d9-b4e9-4eca-b847-748d6c7bafa9"
 Metaheuristics = "bcdb8e00-2c21-11e9-3065-2b553b22f898"
 MillerExtendedHarmonic = "c82744c2-dc08-461a-8c37-87ab04d0f9b8"
+NEO = "081d0416-0f87-42c0-8034-5895805e76d8"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 NNeutronics = "a9424c20-d414-11ec-167b-9106c24d956c"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"

src/actors/transport/flux_calculator_actor.jl

@@ -36,17 +36,13 @@ function ActorFluxCalculator(dd::IMAS.dd, par::FUSEparameters__ActorFluxCalculat
     elseif par.turbulence_model == :TGLF
         act.ActorTGLF.rho_transport = par.rho_transport
         turb_actor = ActorTGLF(dd, act.ActorTGLF)
-    else
-        error("turbulence_model `$(par.turbulence_model)` is not supported yet")
     end
 
     if par.neoclassical_model == :none
         logging(Logging.Debug, :actors, "ActorFluxCalculator: neoclassical transport disabled")
     elseif par.neoclassical_model == :neoclassical
         act.ActorNeoclassical.rho_transport = par.rho_transport
         neoc_actor = ActorNeoclassical(dd, act.ActorNeoclassical)
-    else
-        error("neoclassical_model `$(par.neoclassical_model)` is not supported yet")
     end
 
     return ActorFluxCalculator(dd, par, turb_actor, neoc_actor)

src/actors/transport/neoclassical_actor.jl

@@ -1,19 +1,19 @@
-import TGLFNN
 import TAUENN
-
+import NEO
 #= ===================== =#
 #  ActorNeoclassical      #
 #= ===================== =#
 Base.@kwdef mutable struct FUSEparameters__ActorNeoclassical{T} <: ParametersActor where {T<:Real}
     _parent::WeakRef = WeakRef(nothing)
     _name::Symbol = :not_set
-    neoclassical_model::Switch{Symbol} = Switch{Symbol}([:changhinton], "-", "Neoclassical model to run"; default=:changhinton)
+    model::Switch{Symbol} = Switch{Symbol}([:changhinton, :neo], "-", "Neoclassical model to run"; default=:changhinton)
     rho_transport::Entry{AbstractVector{<:T}} = Entry{AbstractVector{<:T}}("-", "rho_tor_norm values to compute neoclassical fluxes on"; default=0.2:0.1:0.8)
 end
 
 mutable struct ActorNeoclassical{D,P} <: PlasmaAbstractActor
     dd::IMAS.dd{D}
     par::FUSEparameters__ActorNeoclassical{P}
+    input_neos::Union{Vector{<:NEO.InputNEO},Missing}
     flux_solutions::Vector{<:IMAS.flux_solution}
 end
 
@@ -31,7 +31,12 @@ end
 
 function ActorNeoclassical(dd::IMAS.dd, par::FUSEparameters__ActorNeoclassical; kw...)
     par = par(kw...)
-    return ActorNeoclassical(dd, par, IMAS.flux_solution[])
+    if par.model == :neo
+        input_neos = Vector{NEO.InputNEO}(undef, length(par.rho_transport))
+    else
+        input_neos = missing
+    end
+    return ActorNeoclassical(dd, par, input_neos, IMAS.flux_solution[])
 end
 
 """
@@ -45,14 +50,22 @@ function _step(actor::ActorNeoclassical)
     model = resize!(dd.core_transport.model, :neoclassical; wipe=false)
     m1d = resize!(model.profiles_1d)
     m1d.grid_flux.rho_tor_norm = par.rho_transport
+    cp1d = dd.core_profiles.profiles_1d[]
 
-    if par.neoclassical_model == :changhinton
+    if par.model == :changhinton
         model.identifier.name = "Chang-Hinton"
         eqt = dd.equilibrium.time_slice[]
-        cp1d = dd.core_profiles.profiles_1d[]
-        actor.flux_solutions = [TAUENN.neoclassical_changhinton(eqt, cp1d, rho, 1) for rho in par.rho_transport]
-    else
-        error("$(par.neoclassical_model) is not implemented")
+        actor.flux_solutions = [NEO.changhinton(eqt, cp1d, rho, 1) for rho in par.rho_transport]
+    elseif par.model == :neo 
+        model.identifier.name = "NEO"
+        rho_cp = cp1d.grid.rho_tor_norm
+        gridpoint_cp = [argmin(abs.(rho_cp .- rho)) for rho in par.rho_transport]
+
+        for (idx,i) in enumerate(gridpoint_cp)
+                actor.input_neos[idx] = NEO.InputNEO(dd, i)
+        end
+
+        actor.flux_solutions = [NEO.run_neo(actor.input_neos[idx]) for idx in 1:length(par.rho_transport)]
     end
 
     return actor
@@ -64,17 +77,30 @@ end
 Writes ActorNeoclassical results to dd.core_transport
 """
 function _finalize(actor::ActorNeoclassical)
+    par = actor.par
     dd = actor.dd
     cp1d = dd.core_profiles.profiles_1d[]
     eqt = dd.equilibrium.time_slice[]
 
     model = findfirst(:neoclassical, actor.dd.core_transport.model)
     m1d = model.profiles_1d[]
-    m1d.total_ion_energy.flux = zeros(length(actor.par.rho_transport))
-    for (neoclassical_idx, rho) in enumerate(actor.par.rho_transport)
+
+    if par.model == :neo
+        m1d.electrons.particles.flux = zeros(length(par.rho_transport))
+        m1d.electrons.energy.flux = zeros(length(par.rho_transport))
+    end
+
+    m1d.total_ion_energy.flux = zeros(length(par.rho_transport))
+
+    for (neoclassical_idx, rho) in enumerate(par.rho_transport)
         rho_transp_idx = findfirst(i -> i == rho, m1d.grid_flux.rho_tor_norm)
         rho_cp_idx = argmin(abs.(cp1d.grid.rho_tor_norm .- rho))
         m1d.total_ion_energy.flux[rho_transp_idx] = actor.flux_solutions[neoclassical_idx].ENERGY_FLUX_i * IMAS.gyrobohm_energy_flux(cp1d, eqt)[rho_cp_idx] # W / m^2
+
+        if par.model == :neo
+            m1d.electrons.particles.flux[rho_transp_idx] = actor.flux_solutions[neoclassical_idx].PARTICLE_FLUX_e * IMAS.gyrobohm_particle_flux(cp1d, eqt)[rho_cp_idx]
+            m1d.electrons.energy.flux[rho_transp_idx] = actor.flux_solutions[neoclassical_idx].ENERGY_FLUX_e * IMAS.gyrobohm_energy_flux(cp1d, eqt)[rho_cp_idx]
+        end
     end
 
     return actor