|
| 1 | +module App |
| 2 | + |
| 3 | + using GenieFramework |
| 4 | + using JutulDarcy |
| 5 | + using Jutul |
| 6 | + using HYPRE |
| 7 | + using DataFrames |
| 8 | + using Base64 |
| 9 | + import CairoMakie |
| 10 | + |
| 11 | + @genietools |
| 12 | + |
| 13 | + DEFAULT_CYCLES = 3 |
| 14 | + DEFAULT_DELTA_HOT = 50.0 |
| 15 | + DEFAULT_DELTA_COLD = 30.0 |
| 16 | + DEFAULT_BASE_TEMP = 40.0 |
| 17 | + DEFAULT_CHARGE_PERIOD = 2:5 |
| 18 | + DEFAULT_DISCHARGE_PERIOD = 8:12 |
| 19 | + |
| 20 | + function simulate_hates(; |
| 21 | + t_res = DEFAULT_BASE_TEMP, |
| 22 | + delta_hot = DEFAULT_DELTA_HOT, |
| 23 | + delta_cold = DEFAULT_DELTA_COLD, |
| 24 | + charge_period = DEFAULT_CHARGE_PERIOD, |
| 25 | + discharge_period = DEFAULT_DISCHARGE_PERIOD, |
| 26 | + ncycles = DEFAULT_CYCLES |
| 27 | + ) |
| 28 | + darcy, litre, year, second = si_units(:darcy, :litre, :year, :second) |
| 29 | + |
| 30 | + nx = 100 |
| 31 | + nz = 100 |
| 32 | + |
| 33 | + nx = 20 |
| 34 | + ny = 20 |
| 35 | + |
| 36 | + temperature_top = convert_to_si(40.0, :Celsius) |
| 37 | + pressure_top = convert_to_si(120.0, :bar) |
| 38 | + delta_charge = delta_hot |
| 39 | + delta_discharge = delta_cold |
| 40 | + |
| 41 | + grad_p = 1000*9.81 |
| 42 | + grad_T = 0.3 |
| 43 | + |
| 44 | + # ## Set up the reservoir |
| 45 | + g = CartesianMesh((nx, 1, nz), (250.0, 250.0, 75.0)) |
| 46 | + reservoir = reservoir_domain(g, |
| 47 | + permeability = [0.3, 0.3, 0.1].*darcy, |
| 48 | + porosity = 0.3, |
| 49 | + rock_thermal_conductivity = 2.0, |
| 50 | + fluid_thermal_conductivity = 0.6 |
| 51 | + ) |
| 52 | + |
| 53 | + depth = reservoir[:cell_centroids][3, :] |
| 54 | + # ## Define wells and model |
| 55 | + di = Int(ceil(nx/4)) |
| 56 | + k = Int(ceil(nz/2)) |
| 57 | + Whot = setup_vertical_well(reservoir, 0+di , 1, toe = k, name = :Hot) |
| 58 | + Wcold = setup_vertical_well(reservoir, nx-di+1, 1, toe = k, name = :Cold) |
| 59 | + |
| 60 | + model, parameters = setup_reservoir_model(reservoir, :geothermal, wells = [Whot, Wcold]); |
| 61 | + # ## Set up boundary and initial conditions |
| 62 | + bcells = Int[] |
| 63 | + pressure_res = Float64[] |
| 64 | + temperature_res = Float64[] |
| 65 | + for cell in 1:number_of_cells(g) |
| 66 | + d = depth[cell] |
| 67 | + push!(pressure_res, pressure_top + grad_p*d) |
| 68 | + push!(temperature_res, temperature_top + grad_T*d) |
| 69 | + |
| 70 | + I, J, K = cell_ijk(g, cell) |
| 71 | + if I == 1 || I == nx |
| 72 | + push!(bcells, cell) |
| 73 | + end |
| 74 | + end |
| 75 | + |
| 76 | + bc = flow_boundary_condition(bcells, reservoir, pressure_res[bcells], temperature_res[bcells]) |
| 77 | + # ## Set up the schedule |
| 78 | + |
| 79 | + # ### Set up forces |
| 80 | + |
| 81 | + # We assume we have a supply amounting to 90 C at 25 l/s for storage. During the |
| 82 | + # rest period, we assume the same discharge rate and a temperature of 10 C. |
| 83 | + charge_rate = 25litre/second |
| 84 | + discharge_rate = charge_rate |
| 85 | + temperature_charge = temperature_top + delta_charge |
| 86 | + temperature_discharge = temperature_top - delta_discharge |
| 87 | + |
| 88 | + # Set up forces for charging |
| 89 | + rate_target = TotalRateTarget(charge_rate) |
| 90 | + ctrl_hot = InjectorControl(rate_target, [1.0], density = 1000.0, temperature = temperature_charge) |
| 91 | + rate_target = TotalRateTarget(-charge_rate) |
| 92 | + ctrl_cold = ProducerControl(rate_target) |
| 93 | + forces_charge = setup_reservoir_forces(model, control = Dict(:Hot => ctrl_hot, :Cold => ctrl_cold), bc = bc) |
| 94 | + |
| 95 | + # Set up forces for discharging |
| 96 | + rate_target = TotalRateTarget(discharge_rate) |
| 97 | + ctrl_cold = InjectorControl(rate_target, [1.0], density = 1000.0, temperature = temperature_discharge) |
| 98 | + rate_target = TotalRateTarget(-discharge_rate) |
| 99 | + ctrl_hot = ProducerControl(rate_target) |
| 100 | + forces_discharge = setup_reservoir_forces(model, control = Dict(:Hot => ctrl_hot, :Cold => ctrl_cold), bc = bc) |
| 101 | + |
| 102 | + # ### Set up forces for rest period |
| 103 | + forces_rest = setup_reservoir_forces(model, bc = bc) |
| 104 | + |
| 105 | + # ### Set up timesteps and assign forces to each timestep |
| 106 | + dt = Float64[] |
| 107 | + forces = [] |
| 108 | + month = year/12 |
| 109 | + num_charge = 0 |
| 110 | + num_discharge = 0 |
| 111 | + num_rest = 0 |
| 112 | + for year in 1:ncycles |
| 113 | + for mno in vcat(6:12, 1:5) |
| 114 | + if mno in charge_period |
| 115 | + push!(dt, month) |
| 116 | + push!(forces, forces_charge) |
| 117 | + num_charge += 1 |
| 118 | + elseif mno in discharge_period |
| 119 | + push!(dt, month) |
| 120 | + push!(forces, forces_discharge) |
| 121 | + num_discharge += 1 |
| 122 | + else |
| 123 | + push!(dt, month) |
| 124 | + push!(forces, forces_rest) |
| 125 | + num_rest += 1 |
| 126 | + end |
| 127 | + end |
| 128 | + end |
| 129 | + @info "Set up schedule" num_charge num_discharge num_rest |
| 130 | + # ## Set up initial state |
| 131 | + state0 = setup_reservoir_state(model, Pressure = pressure_res, Temperature = temperature_res) |
| 132 | + # ## Simulate the case |
| 133 | + ws, states, time_s = simulate_reservoir(state0, model, dt, |
| 134 | + forces = forces, |
| 135 | + parameters = parameters, |
| 136 | + info_level = 1 |
| 137 | + ) |
| 138 | + @info "Simulation done" |
| 139 | + # ## Plot energy recovery factor |
| 140 | + # The energy recovery factor η is defined as the amount of stored to produced |
| 141 | + # energy. We plot this both cumulatively and for each of the 25 yearly cycles |
| 142 | + wd = ws.wells[:Hot] |
| 143 | + c_p_water = 4.186 # kJ/kgK |
| 144 | + well_temp = wd[:temperature] |
| 145 | + well_temp_cold = ws.wells[:Cold][:temperature] |
| 146 | + |
| 147 | + q = wd[:mass_rate] |
| 148 | + storage = q .> 0 |
| 149 | + q_store = q.*storage |
| 150 | + q_prod = q.*(.!storage) |
| 151 | + stored_energy = well_temp.*q_store.*c_p_water.*dt |
| 152 | + produced_energy = -well_temp.*q_prod.*c_p_water.*dt |
| 153 | + η_cumulative = cumsum(produced_energy)./cumsum(stored_energy) |
| 154 | + t = cumsum(dt)./si_unit(:day) |
| 155 | + |
| 156 | + @info "Cycles" |
| 157 | + num_years = ncycles |
| 158 | + eta, T = zeros(num_years), zeros(num_years) |
| 159 | + for i = 1:num_years |
| 160 | + ix = (1:12) .+ 12*(i-1) |
| 161 | + se = sum(stored_energy[ix]) |
| 162 | + pe = sum(produced_energy[ix]) |
| 163 | + eta[i] = pe/se |
| 164 | + T[i] = t[ix[end]] |
| 165 | + end |
| 166 | + @info "Cycles done." |
| 167 | + |
| 168 | + return Dict( |
| 169 | + :wtemp => well_temp .- 273.15, |
| 170 | + :ctemp => well_temp_cold .- 273.15, |
| 171 | + :time => time_s./si_unit(:day), |
| 172 | + :T => T, |
| 173 | + :eta => eta, |
| 174 | + :T_spatial => map(x -> reshape(x[:Temperature] .- 273.15, nx, nz), states) |
| 175 | + ) |
| 176 | + end |
| 177 | + |
| 178 | + function figure_to_html(fig) |
| 179 | + # buffer = Base.IOBuffer() |
| 180 | + fname = "tmpfile.png" |
| 181 | + CairoMakie.save(fname, fig) |
| 182 | + buffer = open(fname, "r") |
| 183 | + data = base64encode(buffer) |
| 184 | + close(buffer) |
| 185 | + rm(fname) |
| 186 | + # return html("""<img src="data:image/png;base64,$(data)">""") |
| 187 | + return "data:image/png;base64,$(data)" |
| 188 | + end |
| 189 | + |
| 190 | + @app begin |
| 191 | + @in ncycles = DEFAULT_CYCLES |
| 192 | + @in delta_hot = DEFAULT_DELTA_HOT |
| 193 | + @in delta_cold = DEFAULT_DELTA_COLD |
| 194 | + @in base_temp = DEFAULT_BASE_TEMP |
| 195 | + @in name = "Genie" |
| 196 | + @in start = false |
| 197 | + @in running = false |
| 198 | + @in ButtonProgress_process = false |
| 199 | + @in ButtonProgress_progress = 0.0 |
| 200 | + @in ChargePeriod = RangeData(DEFAULT_CHARGE_PERIOD) |
| 201 | + @in DisChargePeriod = RangeData(DEFAULT_DISCHARGE_PERIOD) |
| 202 | + @in tab_selected = "hot_temp" |
| 203 | + @in tplot_stepno = 0.5 |
| 204 | + @out hotplot = PlotData() |
| 205 | + @out coldplot = PlotData() |
| 206 | + @out etaplot = PlotData() |
| 207 | + @out sim_result = simulate_hates() |
| 208 | + @out imgstr = "" |
| 209 | + @private u_x = [] |
| 210 | + @private u_y = [] |
| 211 | + @onchange tplot_stepno begin |
| 212 | + temperature_spatial = sim_result[:T_spatial] |
| 213 | + nstep = length(temperature_spatial) |
| 214 | + ix = clamp(Int(round(tplot_stepno*nstep)), 1, nstep) |
| 215 | + data = temperature_spatial[ix] |
| 216 | + data = data[:, end:-1:1] |
| 217 | + fig = CairoMakie.Figure(size = (800, 300)) |
| 218 | + ax = CairoMakie.Axis(fig[1, 1], title = "Step $ix/$nstep") |
| 219 | + plt = CairoMakie.heatmap!(ax, data, colormap = :hot, colorrange = (base_temp - delta_cold, base_temp + delta_hot)) |
| 220 | + CairoMakie.Colorbar(fig[1, 2], plt, label = "Temperature / °C") |
| 221 | + imgstr = figure_to_html(fig) |
| 222 | + end |
| 223 | + @onchange ChargePeriod begin |
| 224 | + # Should truncate the discharge period here. |
| 225 | + end |
| 226 | + @onbutton ButtonProgress_process begin |
| 227 | + @info "Hello button clicked" running |
| 228 | + running = false |
| 229 | + # u_x = [] |
| 230 | + # u_y = [] |
| 231 | + empty!(u_x) |
| 232 | + empty!(u_y) |
| 233 | + t = 0.0 |
| 234 | + if running == false |
| 235 | + running = true |
| 236 | + res = simulate_hates( |
| 237 | + t_res = base_temp, |
| 238 | + delta_hot = delta_hot, |
| 239 | + delta_cold = delta_cold, |
| 240 | + charge_period = ChargePeriod.range, |
| 241 | + discharge_period = DisChargePeriod.range, |
| 242 | + ncycles = ncycles |
| 243 | + ) |
| 244 | + @info "Simulation ok? " |
| 245 | + ButtonProgress_progress = 0.0 |
| 246 | + wtemp = res[:wtemp] |
| 247 | + t = res[:time] |
| 248 | + @info "" res |
| 249 | + |
| 250 | + for i in eachindex(wtemp, t) |
| 251 | + push!(u_x, t[i]) |
| 252 | + push!(u_y, wtemp[i]) |
| 253 | + end |
| 254 | + @show ButtonProgress_progress |
| 255 | + hotplot = PlotData( |
| 256 | + x = u_x, |
| 257 | + y = u_y, |
| 258 | + plot = StipplePlotly.Charts.PLOT_TYPE_LINE |
| 259 | + ) |
| 260 | + coldplot = PlotData( |
| 261 | + x = u_x, |
| 262 | + y = res[:ctemp], |
| 263 | + plot = StipplePlotly.Charts.PLOT_TYPE_LINE |
| 264 | + ) |
| 265 | + yr = res[:T] |
| 266 | + eta = res[:eta] |
| 267 | + etaplot = PlotData( |
| 268 | + x = yr, |
| 269 | + y = eta, |
| 270 | + plot = StipplePlotly.Charts.PLOT_TYPE_LINE |
| 271 | + ) |
| 272 | + end |
| 273 | + end |
| 274 | + end |
| 275 | + |
| 276 | + |
| 277 | + function ui() |
| 278 | + [ |
| 279 | + h1("High-temperature aquifer thermal energy storage (HT-ATES)") |
| 280 | + p("Fast simulation of energy storage with Fimbul+JutulDarcy.jl - View the results in the tabs below") |
| 281 | + [ |
| 282 | + tabgroup( |
| 283 | + :tab_selected, |
| 284 | + inlinelabel = true, |
| 285 | + class = "bg-primary text-white shadow-2", |
| 286 | + [ |
| 287 | + tab(name = "hot_temp", icon = "local_fire_department", label = "Hot well"), |
| 288 | + tab(name = "cold_temp", icon = "bolt", label = "Cold well"), |
| 289 | + tab(name = "energy", icon = "bolt", label = "Energy recovered"), |
| 290 | + tab(name = "reservoir", icon = "volcano", label = "Reservoir"), |
| 291 | + ], |
| 292 | + ), |
| 293 | + tabpanels( |
| 294 | + :tab_selected, |
| 295 | + animated = true, |
| 296 | + var"transition-prev" = "scale", |
| 297 | + var"transition-next" = "scale", |
| 298 | + [ |
| 299 | + tabpanel(name = "hot_temp", [ |
| 300 | + plot(:hotplot, layout = PlotLayout( |
| 301 | + xaxis = [ |
| 302 | + PlotLayoutAxis(xy="x", title="Time / days") |
| 303 | + ], |
| 304 | + yaxis = [ |
| 305 | + PlotLayoutAxis(xy="y", title="Temperature / °C") |
| 306 | + ], |
| 307 | + ) |
| 308 | + ) |
| 309 | + ]), |
| 310 | + tabpanel(name = "cold_temp", [ |
| 311 | + plot(:coldplot, layout = PlotLayout( |
| 312 | + xaxis = [ |
| 313 | + PlotLayoutAxis(xy="x", title="Time / years") |
| 314 | + ], |
| 315 | + yaxis = [ |
| 316 | + PlotLayoutAxis(xy="y", title="Temperature / °C") |
| 317 | + ], |
| 318 | + ) |
| 319 | + ) |
| 320 | + ]), |
| 321 | + tabpanel(name = "energy", [ |
| 322 | + plot(:etaplot, layout = PlotLayout( |
| 323 | + xaxis = [ |
| 324 | + PlotLayoutAxis(xy="x", title="Time / years") |
| 325 | + ], |
| 326 | + yaxis = [ |
| 327 | + PlotLayoutAxis(xy="y", title="Produced energy / kJ") |
| 328 | + ], |
| 329 | + ) |
| 330 | + ) |
| 331 | + ]), |
| 332 | + tabpanel(name = "reservoir", [ |
| 333 | + p("Reservoir temperature"), |
| 334 | + # html("{{imgstr}}"), |
| 335 | + imageview(src = :imgstr), |
| 336 | + itemsection(slider(0.0:0.001:1.0, :tplot_stepno, label = "", color = "red")), |
| 337 | + ]), |
| 338 | + ], |
| 339 | + ), |
| 340 | + ] |
| 341 | + |
| 342 | + p("Charging and discharging temperature difference") |
| 343 | + item([ |
| 344 | + itemsection(avatar = "", icon("local_fire_department", color = "red")), |
| 345 | + itemsection(slider(0:1:50, :delta_hot, label = "", color = "red")), |
| 346 | + itemsection(avatar = "", icon("ac_unit", color = "blue")), |
| 347 | + itemsection(slider(0:1:50, :delta_cold, label = "", color = "blue")), |
| 348 | + ]) |
| 349 | + p("Number of yearly cycles to simulate") |
| 350 | + item( |
| 351 | + [ |
| 352 | + itemsection(avatar = "", icon("keyboard_double_arrow_up", color = "teal")), |
| 353 | + itemsection(slider(1:1:50, :ncycles, label = "", color = "teal")) |
| 354 | + ] |
| 355 | + ) |
| 356 | + p("Simulating {{ncycles}} cycles, charge ΔT of {{delta_hot}}°C and discharge ΔT of {{delta_cold}}°C", class="st-module") |
| 357 | + p("Charging period") |
| 358 | + item( |
| 359 | + [ |
| 360 | + itemsection(avatar = "", icon("keyboard_double_arrow_down", color = "red")), |
| 361 | + range(1:1:12, :ChargePeriod, markers = true, label = true, color = "red"), |
| 362 | + ] |
| 363 | + ) |
| 364 | + p("Discharge period") |
| 365 | + item( |
| 366 | + [ |
| 367 | + itemsection(avatar = "", icon("keyboard_double_arrow_up", color = "blue")), |
| 368 | + range(1:1:12, :DisChargePeriod, markers = true, label = true, color = "blue") |
| 369 | + ] |
| 370 | + ) |
| 371 | + |
| 372 | + btn( |
| 373 | + "Run simulation", |
| 374 | + @click(:ButtonProgress_process), |
| 375 | + loading = :ButtonProgress_process, |
| 376 | + percentage = :ButtonProgress_progress, |
| 377 | + color = "primary", |
| 378 | + icon = "rocket_launch", |
| 379 | + class = "q-mr-sm", |
| 380 | + ) |
| 381 | + separator() |
| 382 | + ] |
| 383 | + end |
| 384 | + |
| 385 | + @page("/", ui) |
| 386 | +end |
0 commit comments