Added a Julia ODE example

CaseStudies/noPCM/src/Julia/Calculations.jl

@@ -0,0 +1,54 @@
+module Calculations
+
+import DifferentialEquations as D
+
+function func_V_W(V_tank::Float64)
+    return V_tank
+end
+
+function func_m_W(rho_W::Float64, V_W::Float64)
+    return V_W * rho_W
+end
+
+function func_tau_W(C_W::Float64, h_C::Float64, A_C::Float64, m_W::Float64)
+    return m_W * C_W / (h_C * A_C)
+end
+
+""" Calculates temperature of the water: the average kinetic energy of the particles within the water (degreeC)
+    - Parameter T_C: temperature of the heating coil: the average kinetic energy of the particles within the coil (degreeC)
+    - Parameter T_init: initial temperature: the temperature at the beginning of the simulation (degreeC)
+    - Parameter t_final: final time: the amount of time elapsed from the beginning of the simulation to its conclusion (s)
+    - Parameter A_tol: absolute tolerance
+    - Parameter R_tol: relative tolerance
+    - Parameter t_step: time step for simulation: the finite discretization of time used in the numerical method for solving the computational model (s)
+    - Parameter tau_W: ODE parameter for water related to decay time: derived parameter based on rate of change of temperature of water (s)
+    - Returns: temperature of the water: the average kinetic energy of the particles within the water (degreeC)
+"""
+function func_T_W(T_C::Float64, T_init::Float64, T_final::Float64, A_tol::Float64, R_tol::Float64, t_step::Float64, tau_W::Float64)
+    """ Calculates the rate of change of T_W based on its current value
+        - Parameter T_W: current value of T_W
+        - Parameter _: params (unused)
+        - Parameter t: time(?) (unused)
+        - Returns: derivative of T_W
+    """
+    function f(T_W, _, _)
+        return [-(1.0 / tau_W) * T_W[1] + 1.0 / tau_W * T_C]
+    end
+
+    # Need to specify total time range.
+    t_span = (0.0, T_final)
+    # Need to give derivative function, initial value, and total time range.
+    r = D.ODEProblem(f, [T_init], t_span)
+    # Need to give problem, algorithm (if you want to choose it), relative tolerance,
+    # absolute tolerance, and length of time steps.
+    sol = D.solve(r, alg=D.DP5(), reltol=R_tol, abstol=A_tol, saveat=t_step)
+    # Some post-processing to unwrap the data we want.
+    T_W = Array{Float64}([])
+    for x in sol.u
+        append!(T_W, x[1])
+    end
+
+    return T_W
+end
+
+end

CaseStudies/noPCM/src/Julia/Constants.jl

@@ -0,0 +1,17 @@
+module Constants
+
+const PI = 3.14159265
+const L_MIN = 0.1
+const L_MAX = 50.0
+const RHO_W_MIN = 950.0
+const RHO_W_MAX = 1000.0
+const A_C_MAX = 100000.0
+const C_W_MIN = 4170.0
+const C_W_MAX = 4210.0
+const H_C_MIN = 10.0
+const H_C_MAX = 10000.0
+const T_FINAL_MAX = 86400.0
+const AR_MIN = 1.0e-2
+const AR_MAX = 100.0
+
+end

CaseStudies/noPCM/src/Julia/Control.jl

@@ -0,0 +1,20 @@
+module Control
+
+include("Calculations.jl")
+import .Calculations
+include("InputParameters.jl")
+import .InputParameters
+include("OutputFormat.jl")
+import .OutputFormat
+
+filename = ARGS[1]
+A_C, C_W, h_C, T_init, t_final, L, T_C, t_step, rho_W, D, A_tol, R_tol, E_W = InputParameters.get_input(filename)
+V_tank = InputParameters.derived_values(D, L)
+InputParameters.input_constraints(A_C, C_W, h_C, T_init, t_final, L, T_C, t_step, rho_W, D, E_W)
+V_W = Calculations.func_V_W(V_tank)
+m_W = Calculations.func_m_W(rho_W, V_W)
+tau_W = Calculations.func_tau_W(C_W, h_C, A_C, m_W)
+T_W = Calculations.func_T_W(T_C, T_init, t_final, A_tol, R_tol, t_step, tau_W)
+OutputFormat.write_output(E_W, T_W)
+
+end

CaseStudies/noPCM/src/Julia/InputParameters.jl

@@ -0,0 +1,239 @@
+module InputParameters
+
+include("Constants.jl")
+import .Constants
+
+function get_input(filename::String)
+    infile = open(filename, "r")
+    readline(infile)
+    A_C = parse(Float64, readline(infile))
+    readline(infile)
+    C_W = parse(Float64, readline(infile))
+    readline(infile)
+    h_C = parse(Float64, readline(infile))
+    readline(infile)
+    T_init = parse(Float64, readline(infile))
+    readline(infile)
+    t_final = parse(Float64, readline(infile))
+    readline(infile)
+    L = parse(Float64, readline(infile))
+    readline(infile)
+    T_C = parse(Float64, readline(infile))
+    readline(infile)
+    t_step = parse(Float64, readline(infile))
+    readline(infile)
+    rho_W = parse(Float64, readline(infile))
+    readline(infile)
+    D = parse(Float64, readline(infile))
+    readline(infile)
+    A_tol = parse(Float64, readline(infile))
+    readline(infile)
+    R_tol = parse(Float64, readline(infile))
+    readline(infile)
+    E_W = parse(Float64, readline(infile))
+    close(infile)
+
+    return A_C, C_W, h_C, T_init, t_final, L, T_C, t_step, rho_W, D, A_tol, R_tol, E_W
+end
+
+function derived_values(D::Float64, L::Float64)
+    V_tank = Constants.PI * (D / 2.0) ^ 2.0 * L
+
+    return V_tank
+end
+
+function input_constraints(A_C::Float64, C_W::Float64, h_C::Float64, T_init::Float64, t_final::Float64, L::Float64, T_C::Float64, t_step::Float64, rho_W::Float64, D::Float64, E_W::Float64)
+    if !(A_C <= Constants.Constants.A_C_MAX)
+        print("Warning: ")
+        print("A_C has value ")
+        print(A_C)
+        print(", but is suggested to be ")
+        print("below ")
+        print(Constants.Constants.A_C_MAX)
+        print(" (A_C_max)")
+        println(".")
+    end
+    if !(Constants.Constants.C_W_MIN < C_W && C_W < Constants.Constants.C_W_MAX)
+        print("Warning: ")
+        print("C_W has value ")
+        print(C_W)
+        print(", but is suggested to be ")
+        print("between ")
+        print(Constants.Constants.C_W_MIN)
+        print(" (C_W_min)")
+        print(" && ")
+        print(Constants.Constants.C_W_MAX)
+        print(" (C_W_max)")
+        println(".")
+    end
+    if !(Constants.Constants.H_C_MIN <= h_C && h_C <= Constants.Constants.H_C_MAX)
+        print("Warning: ")
+        print("h_C has value ")
+        print(h_C)
+        print(", but is suggested to be ")
+        print("between ")
+        print(Constants.Constants.H_C_MIN)
+        print(" (h_C_min)")
+        print(" && ")
+        print(Constants.Constants.H_C_MAX)
+        print(" (h_C_max)")
+        println(".")
+    end
+    if !(t_final < Constants.Constants.T_FINAL_MAX)
+        print("Warning: ")
+        print("t_final has value ")
+        print(t_final)
+        print(", but is suggested to be ")
+        print("below ")
+        print(Constants.Constants.T_FINAL_MAX)
+        print(" (t_final_max)")
+        println(".")
+    end
+    if !(Constants.Constants.L_MIN <= L && L <= Constants.Constants.L_MAX)
+        print("Warning: ")
+        print("L has value ")
+        print(L)
+        print(", but is suggested to be ")
+        print("between ")
+        print(Constants.Constants.L_MIN)
+        print(" (L_min)")
+        print(" && ")
+        print(Constants.Constants.L_MAX)
+        print(" (L_max)")
+        println(".")
+    end
+    if !(Constants.Constants.RHO_W_MIN < rho_W && rho_W <= Constants.Constants.RHO_W_MAX)
+        print("Warning: ")
+        print("rho_W has value ")
+        print(rho_W)
+        print(", but is suggested to be ")
+        print("between ")
+        print(Constants.Constants.RHO_W_MIN)
+        print(" (rho_W_min)")
+        print(" && ")
+        print(Constants.Constants.RHO_W_MAX)
+        print(" (rho_W_max)")
+        println(".")
+    end
+    if !(Constants.Constants.AR_MIN <= D && D <= Constants.Constants.AR_MAX)
+        print("Warning: ")
+        print("D has value ")
+        print(D)
+        print(", but is suggested to be ")
+        print("between ")
+        print(Constants.Constants.AR_MIN)
+        print(" (AR_min)")
+        print(" && ")
+        print(Constants.Constants.AR_MAX)
+        print(" (AR_max)")
+        println(".")
+    end
+    if !(A_C > 0.0)
+        print("Warning: ")
+        print("A_C has value ")
+        print(A_C)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(C_W > 0.0)
+        print("Warning: ")
+        print("C_W has value ")
+        print(C_W)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(h_C > 0.0)
+        print("Warning: ")
+        print("h_C has value ")
+        print(h_C)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(0.0 < T_init && T_init < 100.0)
+        print("Warning: ")
+        print("T_init has value ")
+        print(T_init)
+        print(", but is suggested to be ")
+        print("between ")
+        print(0.0)
+        print(" && ")
+        print(100.0)
+        println(".")
+    end
+    if !(t_final > 0.0)
+        print("Warning: ")
+        print("t_final has value ")
+        print(t_final)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(L > 0.0)
+        print("Warning: ")
+        print("L has value ")
+        print(L)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(0.0 < T_C && T_C < 100.0)
+        print("Warning: ")
+        print("T_C has value ")
+        print(T_C)
+        print(", but is suggested to be ")
+        print("between ")
+        print(0.0)
+        print(" && ")
+        print(100.0)
+        println(".")
+    end
+    if !(0.0 < t_step && t_step < t_final)
+        print("Warning: ")
+        print("t_step has value ")
+        print(t_step)
+        print(", but is suggested to be ")
+        print("between ")
+        print(0.0)
+        print(" && ")
+        print(t_final)
+        print(" (t_final)")
+        println(".")
+    end
+    if !(rho_W > 0.0)
+        print("Warning: ")
+        print("rho_W has value ")
+        print(rho_W)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(D > 0.0)
+        print("Warning: ")
+        print("D has value ")
+        print(D)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+    if !(E_W >= 0.0)
+        print("Warning: ")
+        print("E_W has value ")
+        print(E_W)
+        print(", but is suggested to be ")
+        print("above ")
+        print(0.0)
+        println(".")
+    end
+end
+
+end

CaseStudies/noPCM/src/Julia/Makefile

@@ -0,0 +1,13 @@
+# Generated by Drasil v0.1-alpha
+
+# Project Name: SWHSNoPCM
+
+# Project Purpose: Investigate the heating of water in a solar water heating
+# tank.
+
+build:
+
+run: build
+	julia Control.jl $(RUNARGS)
+
+.PHONY: build run

CaseStudies/noPCM/src/Julia/OutputFormat.jl

@@ -0,0 +1,12 @@
+module OutputFormat
+
+function write_output(E_W::Float64, T_W::Array{Float64})
+    outputfile = open("output.txt", "w")
+    print(outputfile, "T_W = ")
+    println(outputfile, T_W)
+    print(outputfile, "E_W = ")
+    println(outputfile, E_W)
+    close(outputfile)
+end
+
+end

CaseStudies/noPCM/src/Julia/input.txt

@@ -0,0 +1,26 @@
+# heating coil surface area (m^2)
+0.12
+# specific heat capacity of water (J/(kg degreeC))
+4186.0
+# convective heat transfer coefficient between coil and water (W/(m^2 degreeC))
+1000.0
+# initial temperature (degreeC)
+40.0
+# final time (s)
+50000.0
+# length of tank (m)
+1.5
+# temperature of the heating coil (degreeC)
+50.0
+# time step for simulation (s)
+10.0
+# density of water (kg/m^3)
+1000.0
+# diameter of tank (m)
+0.412
+# absolute tolerance
+1.0e-10
+# relative tolerance
+1.0e-10
+# change in heat energy in the water (J)
+0.0