feedback_control_scenario_workflow

[nanshi1177]

ini,act = FUSE.case_parameters(:FPP)
ini.core_profiles.ne_setting = :greenwald_fraction_ped
ini.core_profiles.ne_value = 0.8
dd = IMAS.dd()
FUSE.init(dd, ini, act; do_plot=false);

act.ActorTEQUILA.relax = 0.1
act.ActorTEQUILA.number_of_iterations=100

ini.pellet_launcher = FUSE.ParametersInits(;n_pl=1).pellet_launcher
ini.pellet_launcher[1].frequency = 8
ini.pellet_launcher[1].size = [2e-3]
ini.pellet_launcher[1].shape = :spherical
ini.pellet_launcher[1].species = :DT

ini.ec_launcher = FUSE.ParametersInits(;n_ec=3).ec_launcher
ini.ec_launcher[1].power_launched = 6.5e+07
ini.ec_launcher[1].rho_0 = 0.65
ini.ec_launcher[1].width = 0.05
ini.ec_launcher[2].power_launched = 0.0e+07
ini.ec_launcher[2].rho_0 = 0.55
ini.ec_launcher[2].width = 0.05
ini.ec_launcher[3].power_launched = 5.0e+05
ini.ec_launcher[3].rho_0 = 0.0
ini.ec_launcher[3].width = 0.05

act.ActorHCD.ec_model = :ECsimple
act.ActorHCD.pellet_model = :Pelletsimple
dd = IMAS.dd()
FUSE.init(dd,ini,act)
chk = FUSE.Checkpoint()
chk[:init_fpp] = dd,ini,act

using NumericalIntegration
my_current = "/home/shinan/.julia/dev/FUSE/playground/shi/highbetap_SAT0_EM_TJLF_fixIp/"
#dd = IMAS.json2imas(my_current * "Iteration_39")
dd_orig = deepcopy(dd)
json_num = 41
betaN_target = 4.2
f_Gr_target = 1.3
diff_max_betaN = 0.05
diff_max_density = 0.05
coef_source_decrease = 0.5
coef_diff_betaN = 0.5
density_err_limit = 2.5
act.ActorTGLF.sat_rule = :sat0
act.ActorTGLF.lump_ions = false
act.ActorFluxMatcher.optimizer_algorithm =:simple
act.ActorFluxMatcher.max_iterations = 40
act.ActorFluxMatcher.rho_transport = 0.2:0.1:0.8
act.ActorFluxMatcher.step_size = 0.1
act.ActorFluxMatcher.evolve_densities=:flux_match
act.ActorTGLF.model = :TJLF
act.ActorNeoclassical.model = :hirshmansigmar
num_iterations = 1
for iteration in 1:num_iterations
    println("========== Start of Iteration $iteration ==========")
    num=0
    diff_betaN = 1.0
    dd0 = deepcopy(dd)
    while abs(diff_betaN) > diff_max_betaN
        FUSE.ActorHCD(dd,act)
        FUSE.ActorFluxMatcher(dd,act;verbose=true,do_plot=true)
        betaN = dd.core_profiles.global_quantities.beta_tor_norm
        diff_betaN = - (betaN[1] - betaN_target) / betaN_target
        Paux_scale = 1.0 + coef_diff_betaN * diff_betaN
        if Paux_scale < 0.9
            Paux_scale = 0.9
        elseif Paux_scale > 1.1
            Paux_scale = 1.1
        end  
        num=num+1
        dd.pulse_schedule.ec.beam[1].power_launched.reference *= Paux_scale
        dd.pulse_schedule.ec.beam[2].power_launched.reference *= Paux_scale  
        println("betaN = ", dd.core_profiles.global_quantities.beta_tor_norm[1])
        println("diff_betaN = ",diff_betaN)
        println("Paux_scale = ",Paux_scale)
        println("Paux_power[1] = ",dd.pulse_schedule.ec.beam[1].power_launched.reference)
        println("Paux_power[2] = ",dd.pulse_schedule.ec.beam[2].power_launched.reference)
        println("Number of Transport Iteration = ",num)
        if num >1 
            println("The number of tranport iteration is beyond the limit.")
            break
        elseif dd.core_profiles.profiles_1d[1].electrons.density[1]/1.e20 > density_err_limit
            dd = deepcopy(dd0)
            dd.pulse_schedule.pellet.launcher[1].frequency.reference *= 0.9
            println("The density ne[0] is beyond the limit.")
            break
        end
        sleep(15)
    end    
    a = 1.375
    Ip = dd.equilibrium.time_slice[].global_quantities.ip / 1e6
    n_Gr = Ip / (pi * a ^2)
    rho_TGYRO_output = dd.core_profiles.profiles_1d[1].grid.rho_tor_norm
    ne_TGYRO_output = dd.core_profiles.profiles_1d[1].electrons.density
    ne_line_avg = integrate(rho_TGYRO_output, ne_TGYRO_output) / 1e20
    f_Gr = ne_line_avg / n_Gr
    density_ratio = f_Gr / f_Gr_target
    ext_particle_coef = 1- (density_ratio-1) * coef_source_decrease
    if abs(density_ratio-1) > diff_max_density
        if ext_particle_coef > 1.1
            ext_particle_coef = 1.1
        elseif ext_particle_coef < 0.9
            ext_particle_coef =0.9
        end
    end
    println("f_Gr = ", f_Gr)
    println("density_ratio = ", density_ratio)
    println("ext_particle_coef = ",ext_particle_coef)
    dd.pulse_schedule.pellet.launcher[1].frequency.reference *= ext_particle_coef 
    println("pellet_launcher.frequency = ",dd.pulse_schedule.pellet.launcher[1].frequency.reference)
    
    #------------ for standard H-mode scenario ------------------
    #FUSE.ActorEquilibrium(dd,act;model=:TEQUILA,ip_from = :pulse_schedule)
    #FUSE.ActorCurrent(dd,act;ip_from = :pulse_schedule)
    #println("betaN = ",dd.core_profiles.global_quantities.beta_tor_norm)
    #println("Ip = ",dd.core_profiles.global_quantities.ip)
    #------------------------------------------------------------
    
    #------------ for high-betap scenario -----------------------
    oh_com = 0.05
    foh = 1.0
    target_foh = 0.2
    relax = 0.2
    while abs(foh - target_foh) > oh_com
        FUSE.ActorCurrent(dd,act;ip_from = :pulse_schedule)
        curoh = dd.core_profiles.global_quantities.ip[1] - dd.core_profiles.global_quantities.current_non_inductive[1]
        foh = curoh / dd.core_profiles.global_quantities.ip[1]
        println("foh = ",foh)
        if abs(foh - target_foh) > oh_com
            dd.pulse_schedule.flux_control.i_plasma.reference *= (1. - relax * (foh - target_foh))
            FUSE.ActorEquilibrium(dd,act;model = :TEQUILA,ip_from = :pulse_schedule)
            println("foh - target_foh = ",foh - target_foh)
            println("i_plasma = ", dd.pulse_schedule.flux_control.i_plasma.reference)
        end
    end
    #--------------------------------------------------------------
    
    json_num += 1
    println("json_num = ",json_num)
    IMAS.imas2json(dd, my_current * "Iteration_$(json_num)";freeze=false)
    println("========== End of Iteration $iteration ==========\n")
    sleep(15)
end

[TimSlendebroek]

We discussed doing this more in the MOOP way by setting constraints and targeting convergence as objective let's work on it monday :)