<commit_before>import "compiler.liszt"
local Grid  = terralib.require 'compiler.grid'
local cmath = terralib.includecstring [[ #include <math.h> ]]


-----------------------------------------------------------------------------
--[[                             OPTIONS                                 ]]--
-----------------------------------------------------------------------------


local grid_options = {
    xnum = 16,
    ynum = 16,
    pos = {0.0, 0.0},
    width = 6.28,
    height = 6.28
}


local particle_options = {
    num = 50
}


-----------------------------------------------------------------------------
--[[                       FLOW/ PARTICLE RELATIONS                      ]]--
-----------------------------------------------------------------------------


-- Declare and initialize grid  and related fields

local grid = Grid.New2dUniformGrid(grid_options.xnum, grid_options.ynum,
                                   grid_options.pos,
                                   grid_options.width, grid_options.height)

-- conserved variables
grid.cells:NewField('rho', L.float)
grid.cells:NewField('rho_velocity', L.vec2f)
grid.cells:NewField('rho_energy', L.float)

-- primitive variables
grid.cells:NewField('velocity', L.vec2f)
grid.cells:NewField('temperature', L.float)
grid.cells:NewField('pressure', L.float)

-- scratch (temporary) fields
-- intermediate value
grid.cells:NewField('rho_temp', L.float)
grid.cells:NewField('rho_velocity_temp', L.vec2f)
grid.cells:NewField('rho_energy_temp', L.float)
-- derivatives
grid.cells:NewField('rho', L.float)
grid.cells:NewField('rho_velocity_t', L.vec2f)
grid.cells:NewGield('rho_energy_t', L.float)
-- flux
-- TODO: Define face and related macros
grid.faces:NewField('rho_flux', L.float)
grid.faces:NewField('rho_velocity_flux', L.vec2f)
grid.faces:NewField('rho_energy_flux', L.float)
grid.faces:NewField('rho_enthalpy', L.float)


-- Declare and initialize particle relation and fields over the particle

local particles = L.NewRelation(particle_options.num, 'particles')

particles:NewField('dual_cell', grid.dual_cells)
particles:NewField('position', L.vec2f)
particles:NewField('velocity', L.vec2f)
particles:NewField('temperature', L.float)

particles:NewField('diameter', L.float)
particles:NewField('density', L.float)

particles:NewField('delta_velocity_over_relaxaion_time', L.vec2f)
particles:NewField('delta_temperature_term', L.float)

-- scratch (temporary) fields
-- derivatives
particles:NewField('position_t', L.vec2f)
particles:NewField('velocity_t', L.vec2f)
particles:NewField('temperature_t', L.float)


-----------------------------------------------------------------------------
--[[                            LISZT GLOBALS                            ]]--
-----------------------------------------------------------------------------


-- If you want to use global variables like scalar values or vectors within
-- a Liszt kernel (variables that are not constants), declare them here.
-- If you use lua variables, Liszt will treat them as constant values.


-----------------------------------------------------------------------------
--[[                             LISZT MACROS                            ]]--
-----------------------------------------------------------------------------


local Rho = L.NewMacro(function(c)
    return liszt `c.rho
end)

local Velocity = L.NewMacro(function(c)
    return liszt `c.velocity
end)

local Temperature = L.NewMacro(function(c)
    return liszt `c.temperature
end)

local InterpolateBilinear = L.NewMacro(function(dc, Field)
    return liszt quote
        var c00 = dc.cell(-1, -1)
        var c01 = dc.cell(-1, 1)
        var c10 = dc.cell(1, -1)
        var c11 = dc.cell(1, 1)
        var delta1 = cmath.abs(c00.position[1] - dc.y)
        var delta2 = cmath.abs(c11.position[1] - dc.y)
        var f1 = (delta1*Field(c00) + delta2*Field(c01)) / (delta1 + delta2)
        var f2 = (delta1*Field(c10) + delta2*Field(c11)) / (delta1 + delta2)
        delta1 = cmath.abs(c00.position[0] - dc.x)
        delta2 = cmath.abs(c11.position[0] - dc.x)
    in
        (delta1*f1 + delta2*f2) / (delta1 + delta2)
    end
end)

local GetDynamicViscosity = L.NewMacro(function(temperature)
    return liszt `fluid_options.dynamic_viscosity_ref *
        cmath.pow(temperature/fluid_options.dynamic_viscosity_temp_ref, 0.75)
end)

-----------------------------------------------------------------------------
--[[                            LISZT KERNELS                            ]]--
-----------------------------------------------------------------------------


-- Locate particles
local LocateParticles = liszt kernel(p : particles)
    p.dual_cell = grid.dual(p.position)
end

-- Initialize enthalpy
local AddInviscidInitialize = liszt kernel(c : grid.cells)
    c.rho_enthalpy = c.rho_energy + c.pressure
end


-- Interpolation for flux values
-- read conserved variables, write flux variables
local AddInviscidGetFlux = liszt kernel(f : grid.faces)
    var num_coeffs = spatial_stencil.num_interpolate_coeffs
    var coeffs     = spatial_stencil.interpolate_coeffs
    var rho_diagonal = 0.0
    var rho_skew     = 0.0
    var rho_velocity_diagonal = {0.0, 0.0}
    var rho_velocity_skew     = {0.0, 0.0}
    var rho_energy_diagonal   = 0.0
    var rho_energy_skew       = 0.0
    var fpdiag = 0.0
    var axis = f.axis
    for i = 0, num_coeffs do
        rho_diagonal += coeffs[i] *
                      ( f.cell(i+1).rho *
                        f.cell(i+1).velocity[axis] +
                        f.cell(i-1).rho *
                        f.cell(i-1).velocity[axis] )
        rho_velocity_diagonal += coeffs[i] *
                               ( f.cell(i+1).rho_velocity *
                                 f.cell(i+1).velocity[axis] +
                                 f.cell(i-1).rho_velocity *
                                 f.cell(i-1).velocity[axis] )
        rho_energy_diagonal += coeffs[i] *
                             ( f.cell(i+1).rho_enthalpy *
                               f.cell(i+1).velocity[axis] +
                               f.cell(i-1).rho_enthalpy *
                               f.cell(i-1).velocity[axis] )
        fpdiag += coeffs[i] *
                ( f.cell(i+1).pressure +
                  f.cell(i-1).pressure )
    end
    -- TODO: I don't understand what skew is. Ask Ivan.
    var s = spatial_stencil.split
    f.rho_flux          = s * rho_diagonal +
                          (1-s) * rho_skew
    f.rho_velocity_flux = s * rho_velocity_diagonal +
                          (1-s) * rho_velocity_skew
    f.rho_energy_flux   = s * rho_energy_diagonal +
                          (1-s) * rho_energy_skew
end


-- Update conserved variables using flux values from part 1
-- write conserved variables, read flux variables
local AddInviscidUpdateUsingFlux = liszt kernel(c : grid.cells)
    c.rho_t -= (c.face(1,0,0).rho_flux -
                c.face(-1,0,0).rho_flux)/c.dx
    c.rho_t -= (c.face(0,1,0).rho_flux -
                c.face(1,-1,0.rho_flux))/c.dy
    c.rho_t -= (c.face(0,0,1).rho_flux -
                c.face(0,0,-1).rho_flux)/c.dz
    c.rho_velocity_t -= (c.face(1,0,0).rho_velocity_flux -
                         c.face(-1,0,0).rho_velocity_flux)/c.dx
    c.rho_velocity_t -= (c.face(0,1,0).rho_velocity_flux -
                         c.face(1,-1,0).rho_velocity_flux)/c.dy
    c.rho_velocity_t -= (c.face(0,0,1).rho_velocity_flux -
                         c.face(0,0,-1).rho_velocity_flux)/c.dz
    c.rho_energy_t -= (c.face(1,0,0).rho_energy_flux -
                       c.face(-1,0,0).rho_energy_flux)/c.dx
    c.rho_energy_t -= (c.face(0,1,0).rho_energy_flux -
                       c.face(1,-1,0).rho_energy_flux)/c.dy
    c.rho_energy_t -= (c.face(0,0,1).rho_energy_flux -
                       c.face(0,0,-1).rho_energy_flux)/c.dz
end


-- Update particle fields based on flow fields
local AddFlowCouplingPartOne = liszt kernel(p: particles)
    var dc = p.dual_cell
    var flow_density     = InterpolateBilinear(dc, Rho)
    var flow_velocity    = InterpolateBilinear(dc, Velocity)
    var flow_temperature = InterpolateBilineat(dc, Temperature)
    var flow_dyn_viscosity = GetDynamicViscosity(flow_temperature)
    p.position_t    += p.velocity
    var relaxation_time = p.density * cmath.pow(p.diameter, 2) /
        (18.0 * flow_dyn_viscosity)
    p.delta_velocity_over_relaxation_time = (flow_velocity - p.velocity)/
        relaxation_time
    var particle_mass = 3.14 * p.density * cmath.pow(p.diameter, 3) / 6.0
    p.delta_temperature_term = 3.14 * cmath.pow(p.diameter, 2) *
        particle_options.convective_coefficient *
        (flow_temperature - p.temperature)
end


local AddFlowCouplingPartTwo = liszt kernel(p : particles)
    p.velocity_t += p.delta_velocity_over_relaxation_time
    p.temperature_t += p.delta_temperature_term/
        (particle_mass * particle_options.heat_capacity)
end


-- Update flow variables using derivatives
-- (intermediate steps in a single time step)
local UpdateFlowFieldsIntermediate = liszt kernel(c : grid.cells)
end


-- Update flow variables using derivatives
-- (last step in a single time step)
local UpdateFlowFieldsFinal = liszt kernel(c : grid.cells)
end


-- Update particle variables using derivatives
-- (intermediate steps in a single time step)
local UpdateParticleFieldsIntermediate = liszt kernel(p : particles)
end


-- Update flow variables using derivatives
-- (last step in a single time step)
local UpdateParticleFieldsFinal = liszt kernel(p : particles)
end


local UpdateAuxiliaryGridVelocity = liszt kernel(c : grid.cells)
    c.velocity = c.rho_velocity / c.rho
end


local UpdateAuxiliaryGridThermodynamics = liszt kernel(c : grid.cells)
end


local UpdateAuxiliaryParticles = liszt kernel(p : particles)
end


-----------------------------------------------------------------------------
--[[                             MAIN LOOP                               ]]--
-----------------------------------------------------------------------------

<commit_msg>adding main loop and time integration for PSAAP2 example<commit_after>import "compiler.liszt"
local Grid  = terralib.require 'compiler.grid'
local cmath = terralib.includecstring [[ #include <math.h> ]]


-----------------------------------------------------------------------------
--[[                             OPTIONS                                 ]]--
-----------------------------------------------------------------------------


local grid_options = {
    xnum = 16,
    ynum = 16,
    pos = {0.0, 0.0},
    width = 6.28,
    height = 6.28
}


local particle_options = {
    num = 50
}


local spatial_stencil = {
}


local time_integrator = {
}


-----------------------------------------------------------------------------
--[[                       FLOW/ PARTICLE RELATIONS                      ]]--
-----------------------------------------------------------------------------


-- Declare and initialize grid  and related fields

local grid = Grid.New2dUniformGrid(grid_options.xnum, grid_options.ynum,
                                   grid_options.pos,
                                   grid_options.width, grid_options.height)

-- conserved variables
grid.cells:NewField('rho', L.float)
grid.cells:NewField('rho_velocity', L.vec2f)
grid.cells:NewField('rho_energy', L.float)

-- primitive variables
grid.cells:NewField('velocity', L.vec2f)
grid.cells:NewField('temperature', L.float)
grid.cells:NewField('pressure', L.float)

-- scratch (temporary) fields
-- intermediate value and copies
grid.cells:NewField('rho_copy', L.float)
grid.cells:NewField('rho_velocity_copy', L.vec2f)
grid.cells:NewField('rho_energy_copy', L.float)
grid.cells:NewField('rho_temp', L.float)
grid.cells:NewField('rho_velocity_temp', L.vec2f)
grid.cells:NewField('rho_energy_temp', L.float)
-- derivatives
grid.cells:NewField('rho_t', L.float)
grid.cells:NewField('rho_velocity_t', L.vec2f)
grid.cells:NewGield('rho_energy_t', L.float)
-- flux
-- TODO: Define edge and related macros
grid.edges:NewField('rho_flux', L.float)
grid.edges:NewField('rho_velocity_flux', L.vec2f)
grid.edges:NewField('rho_energy_flux', L.float)
grid.edges:NewField('rho_enthalpy', L.float)


-- Declare and initialize particle relation and fields over the particle

local particles = L.NewRelation(particle_options.num, 'particles')

particles:NewField('dual_cell', grid.dual_cells)
particles:NewField('position', L.vec2f)
particles:NewField('velocity', L.vec2f)
particles:NewField('temperature', L.float)

particles:NewField('diameter', L.float)
particles:NewField('density', L.float)

particles:NewField('delta_velocity_over_relaxaion_time', L.vec2f)
particles:NewField('delta_temperature_term', L.float)

-- scratch (temporary) fields
-- intermediate values and copies
particles:NewField('position_copy', L.vec2f)
particles:NewField('velocity_copy', L.vec2f)
particles:NewField('temperature_copy', L.float)
particles:NewField('position_temp', L.vec2f)
particles:NewField('velocity_temp', L.vec2f)
particles:NewField('temperature_temp', L.float)
-- derivatives
particles:NewField('position_t', L.vec2f)
particles:NewField('velocity_t', L.vec2f)
particles:NewField('temperature_t', L.float)


-----------------------------------------------------------------------------
--[[                            LISZT GLOBALS                            ]]--
-----------------------------------------------------------------------------


-- If you want to use global variables like scalar values or vectors within
-- a Liszt kernel (variables that are not constants), declare them here.
-- If you use lua variables, Liszt will treat them as constant values.

local delta_time = L.NewGlobal(L.float, 0.05)


-----------------------------------------------------------------------------
--[[                             LISZT MACROS                            ]]--
-----------------------------------------------------------------------------


local Rho = L.NewMacro(function(r)
    return liszt `r.rho
end)

local Velocity = L.NewMacro(function(r)
    return liszt `r.velocity
end)

local Temperature = L.NewMacro(function(r)
    return liszt `r.temperature
end)

local InterpolateBilinear = L.NewMacro(function(dc, Field)
    return liszt quote
        var c00 = dc.cell(-1, -1)
        var c01 = dc.cell(-1, 1)
        var c10 = dc.cell(1, -1)
        var c11 = dc.cell(1, 1)
        var delta1 = cmath.abs(c00.position[1] - dc.y)
        var delta2 = cmath.abs(c11.position[1] - dc.y)
        var f1 = (delta1*Field(c00) + delta2*Field(c01)) / (delta1 + delta2)
        var f2 = (delta1*Field(c10) + delta2*Field(c11)) / (delta1 + delta2)
        delta1 = cmath.abs(c00.position[0] - dc.x)
        delta2 = cmath.abs(c11.position[0] - dc.x)
    in
        (delta1*f1 + delta2*f2) / (delta1 + delta2)
    end
end)

local GetDynamicViscosity = L.NewMacro(function(temperature)
    return liszt `fluid_options.dynamic_viscosity_ref *
        cmath.pow(temperature/fluid_options.dynamic_viscosity_temp_ref, 0.75)
end)


-----------------------------------------------------------------------------
--[[                            LISZT KERNELS                            ]]--
-----------------------------------------------------------------------------


-- Locate particles
local LocateParticles = liszt kernel(p : particles)
    p.dual_cell = grid.dual(p.position)
end

-- Initialize enthalpy and derivatives
local AddInviscidInitialize = liszt kernel(c : grid.cells)
    c.rho_enthalpy = c.rho_energy + c.pressure
end


-- Interpolation for flux values
-- read conserved variables, write flux variables
local AddInviscidGetFlux = liszt kernel(e : grid.edges)
    var num_coeffs = spatial_stencil.num_interpolate_coeffs
    var coeffs     = spatial_stencil.interpolate_coeffs
    var rho_diagonal = 0.0
    var rho_skew     = 0.0
    var rho_velocity_diagonal = {0.0, 0.0}
    var rho_velocity_skew     = {0.0, 0.0}
    var rho_energy_diagonal   = 0.0
    var rho_energy_skew       = 0.0
    var epdiag = 0.0
    var axis = f.axis
    for i = 0, num_coeffs do
        rho_diagonal += coeffs[i] *
                      ( e.cell(i+1).rho *
                        e.cell(i+1).velocity[axis] +
                        e.cell(i-1).rho *
                        e.cell(i-1).velocity[axis] )
        rho_velocity_diagonal += coeffs[i] *
                               ( e.cell(i+1).rho_velocity *
                                 e.cell(i+1).velocity[axis] +
                                 e.cell(i-1).rho_velocity *
                                 e.cell(i-1).velocity[axis] )
        rho_energy_diagonal += coeffs[i] *
                             ( e.cell(i+1).rho_enthalpy *
                               e.cell(i+1).velocity[axis] +
                               e.cell(i-1).rho_enthalpy *
                               e.cell(i-1).velocity[axis] )
        epdiag += coeffs[i] *
                ( e.cell(i+1).pressure +
                  e.cell(i-1).pressure )
    end
    -- TODO: I don't understand what skew is. Ask Ivan.
    var s = spatial_stencil.split
    e.rho_flux          = s * rho_diagonal +
                          (1-s) * rho_skew
    e.rho_velocity_flux = s * rho_velocity_diagonal +
                          (1-s) * rho_velocity_skew
    e.rho_energy_flux   = s * rho_energy_diagonal +
                          (1-s) * rho_energy_skew
end


-- Update conserved variables using flux values from previous part
-- write conserved variables, read flux variables
local AddInviscidUpdateUsingFlux = liszt kernel(c : grid.cells)
    c.rho_t -= (c.edge(1,0).rho_flux -
                c.edge(-1,0).rho_flux)/c.dx
    c.rho_t -= (c.edge(0,1).rho_flux -
                c.edge(1,-1.rho_flux))/c.dy
    c.rho_velocity_t -= (c.edge(1,0).rho_velocity_flux -
                         c.edge(-1,0).rho_velocity_flux)/c.dx
    c.rho_velocity_t -= (c.edge(0,1).rho_velocity_flux -
                         c.edge(1,-1).rho_velocity_flux)/c.dy
    c.rho_energy_t -= (c.edge(1,0).rho_energy_flux -
                       c.edge(-1,0).rho_energy_flux)/c.dx
    c.rho_energy_t -= (c.edge(0,1).rho_energy_flux -
                       c.edge(1,-1).rho_energy_flux)/c.dy
end


-- Update particle fields based on flow fields
local AddFlowCouplingPartOne = liszt kernel(p: particles)
    var dc = p.dual_cell
    var flow_density     = InterpolateBilinear(dc, Rho)
    var flow_velocity    = InterpolateBilinear(dc, Velocity)
    var flow_temperature = InterpolateBilineat(dc, Temperature)
    var flow_dyn_viscosity = GetDynamicViscosity(flow_temperature)
    p.position_t    += p.velocity
    var relaxation_time = p.density * cmath.pow(p.diameter, 2) /
        (18.0 * flow_dyn_viscosity)
    p.delta_velocity_over_relaxation_time = (flow_velocity - p.velocity)/
        relaxation_time
    var particle_mass = 3.14 * p.density * cmath.pow(p.diameter, 3) / 6.0
    p.delta_temperature_term = 3.14 * cmath.pow(p.diameter, 2) *
        particle_options.convective_coefficient *
        (flow_temperature - p.temperature)
end


local AddFlowCouplingPartTwo = liszt kernel(p : particles)
    p.velocity_t += p.delta_velocity_over_relaxation_time
    p.temperature_t += p.delta_temperature_term/
        (particle_mass * particle_options.heat_capacity)
end


-- Kernel generator for update flow variables using derivatives
local UpdateFlowKernels = {}
local function GenerateUpdateFlowKernels(relation, stage)
    local coeff_fun  = time_integrator.coeff_function[stage]
    local coeff_time = time_integrator.coeff_time[stage]
    if stage <= 3 then
        return liszt kernel(r : relation)
            r.rho_temp += coeff_fun * delta_time * r.rho_t
            r.rho       = r.rho_copy +
                          coeff_time * delta_time * r.rho_t
            r.rho_velocity_temp += coeff_fun * delta_time * r.rho_velocity_t
            r.rho_velocity       = r.rho_velocity_copy +
                                   coeff_time * delta_time * r.rho_velocity_t
            r.rho_energy_temp += coeff_fun * delta_time * r.rho_energy_t
            r.rho_energy       = r.rho_energy_copy +
                                 coeff_time * delta_time * r.rho_energy_t
        end
    elseif
        return liszt kernel(r : relation)
            r.rho = r.rho_temp +
                    coeff_fun * delta_time * r.rho_t
            r.rho_velocity = r.rho_velocity_temp +
                             coeff_fun * delta_time * r.rho_velocity_t
            r.rho_energy = r.rho_energy_temp +
                           coeff_fun * delta_time * r.rho_energy_t
        end
    end
end
for i = 1, 4 do
    UpdateFlowKernels[i] = GenerateUpdateFlowKernels(grid.cells, i)
end


-- Update particle variables using derivatives
local UpdateParticleKernels = {}
local function GenerateUpdateParticleKernels(relation, stage)
    local coeff_fun  = time_integrator.coeff_function[stage]
    local coeff_time = time_integrator.coeff_time[stage]
    if stage <= 3 then
        return liszt kernel(r : relation)
            r.position_temp += coeff_fun * delta_time * r.position_t
            r.position       = r.position_copy +
                               coeff_time * delta_time * r.position_t
            r.velocity_temp += coeff_fun * delta_time * r.velocity_t
            r.velocity       = r.velocity_copy +
                               coeff_time * delta_time * r.velocity_t
            r.temperature_temp += coeff_fun * delta_time * r.temperature_t
            r.temperature       = r.temperature_copy +
                                  coeff_time * delta_time * r.temperature_t
        end
    elseif
        return liszt kernel(r : relation)
            r.position = r.position_temp +
                         coeff_fun * delta_time * r.position_t
            r.velocity = r.velocity_temp +
                         coeff_fun * delta_time * r.velocity_t
            r.temperature = r.rho_temperature_temp +
                            coeff_fun * delta_time * r.temperature_t
        end
    end
end
for i = 1, 4 do
    UpdateParticleKernels[i] = GenerateUpdateParticleKernels(particles, i)
end


-- Update time function
local function UpdateTime(stage)
    time_integrator.sim_time = time_integrator.sim_time +
                               time_integrator.coeff_time[stage] *
                               delta_time:value()
end


local UpdateAuxiliaryGridVelocity = liszt kernel(c : grid.cells)
    c.velocity = c.rho_velocity / c.rho
end


local UpdateAuxiliaryGridThermodynamics = liszt kernel(c : grid.cells)
end


local UpdateAuxiliaryParticles = liszt kernel(p : particles)
end


-----------------------------------------------------------------------------
--[[                             MAIN LOOP                               ]]--
-----------------------------------------------------------------------------


local function InitializeTemporaries()
end


local function InitializeDerivatives()
end


local function AddInviscid()
    AddInviscidInitialize(grid.cells)
    AddInviscidGetFlux(grid.edges)
    AddInviscidUpdateUsingFlux(grid.cells)
end


local function AddFlowCoupling()
    AddFlowCouplingPartOne(particles)
    AddFlowCouplingPartTwo(particles)
end

local function UpdateFlow(i)
    UpdateFlowKernels[i](grid.cells)
end


local function UpdateParticles(i)
    UpdateParticleKernels[i](particles)
end


while (time_integrator.sim_time < time_integrator.final_time) do
    InitializeTemporaries()
    for stage = 1, 4 do
        AddInviscid()
        AddFlowCoupling()
        UpdateFlow(i)
        UpdateParticles(i)
    end
end
{'fraction_non_alphanumeric': 0.17828964113005855, 'fraction_numerical': 0.015610418257402222, 'mean_word_length': 1.8008791731020555, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4298768', 'n_tokens_mistral': 7659, 'n_tokens_neox': 7098, 'n_words': 1392}