Can you fix this code: --[[
    Height Field Fluid Simulation in Roblox Studio
    Based on "Height Field Fluids" by Matthias Müller

    This script simulates the shallow water equations over a terrain heightmap.
    The simulation includes terrain interaction, fluid dynamics, numerical stability enhancements, and visualization.
]]

-- User-adjustable parameters
local GRID_WIDTH = 100      -- Number of grid cells in x-direction
local GRID_HEIGHT = 100     -- Number of grid cells in y-direction
local DX = 4                -- Cell size in x (studs)
local DY = 4                -- Cell size in y (studs)
local OFFSET_X = 0          -- World offset in x-direction (adjust to match terrain position)
local OFFSET_Z = 0          -- World offset in z-direction (adjust to match terrain position)
local GRAVITY = 100         -- Acceleration due to gravity (studs/s^2)
local CFL_NUMBER = 0.1      -- Courant number for stability (<= 1)
local WATER_ADDITION_RATE = 100  -- Base rate at which water is added at the source (studs per second)
local WATER_SOURCE_RADIUS = 2    -- Radius of water source area (in cells)

-- Stability and physics parameters
local MAX_WATER_HEIGHT = 100      -- Maximum allowable water height
local MAX_SLOPE = 1               -- Maximum allowable slope between cells
local DAMPING_COEFFICIENT = 0.1   -- Damping coefficient for fluxes
local DIFFUSION_COEFFICIENT = 0.05   -- Base diffusion coefficient for water heights
local VISCOSITY_COEFFICIENT = 0.1    -- Viscosity coefficient
local EPSILON = 1e-6              -- Small epsilon to prevent division by zero
local MIN_WATER_DEPTH = 0.01      -- Minimum water depth to consider a cell wet

-- Data structures
local H_Terrain = {}    -- Terrain height
local H_Water = {}      -- Water depth
local E = {}            -- Water surface elevation (h + z)
local Qx = {}           -- Flux in x-direction
local Qy = {}           -- Flux in y-direction
local SlopeX = {}       -- Terrain slope in x-direction
local SlopeY = {}       -- Terrain slope in y-direction
local WaterParts = {}   -- Visual representation of water

-- Initialize data structures
for i = 0, GRID_WIDTH + 3 do  -- Added buffer cells for boundaries
	H_Terrain[i] = {}
	H_Water[i] = {}
	E[i] = {}
	Qx[i] = {}
	Qy[i] = {}
	SlopeX[i] = {}
	SlopeY[i] = {}
	for j = 0, GRID_HEIGHT + 3 do
		H_Terrain[i][j] = 0
		H_Water[i][j] = 0
		E[i][j] = 0
		Qx[i][j] = 0
		Qy[i][j] = 0
		SlopeX[i][j] = 0
		SlopeY[i][j] = 0
	end
end

-- Function to sample terrain height using raycasts
local function SampleTerrainHeight(i, j)
	-- Convert grid indices to world coordinates
	local worldX = (i - 2) * DX + OFFSET_X  -- Adjusted index for buffer
	local worldZ = (j - 2) * DY + OFFSET_Z  -- Adjusted index for buffer
	local origin = Vector3.new(worldX, 1000, worldZ)   -- Start raycast from above
	local direction = Vector3.new(0, -2000, 0)         -- Raycast downwards

	-- Raycast parameters
	local params = RaycastParams.new()
	params.FilterDescendantsInstances = {workspace.Terrain}  -- Whitelist terrain
	params.FilterType = Enum.RaycastFilterType.Whitelist
	params.IgnoreWater = true  -- Ignore water parts

	local rayResult = workspace:Raycast(origin, direction, params)
	if rayResult then
		local terrainHeight = rayResult.Position.Y
		return terrainHeight
	else
		-- Default terrain height if nothing is hit
		return 0
	end
end

-- Generate the terrain heightmap
function GenerateTerrainHeightmap()
	for i = 0, GRID_WIDTH + 3 do
		for j = 0, GRID_HEIGHT + 3 do
			H_Terrain[i][j] = SampleTerrainHeight(i, j)
			H_Water[i][j] = 0
			E[i][j] = H_Terrain[i][j] + H_Water[i][j]
		end
	end
end

-- Water Source
local WaterSource = {
	x = math.floor(GRID_WIDTH / 2) + 1,
	y = math.floor(GRID_HEIGHT / 2) + 1,
	rate = WATER_ADDITION_RATE  -- Water addition rate
}

function AddWater(deltaTime)
	local radius = WATER_SOURCE_RADIUS
	local totalCells = 0

	-- First, count the total number of cells within the radius
	for i = WaterSource.x - radius, WaterSource.x + radius do
		for j = WaterSource.y - radius, WaterSource.y + radius do
			local dist = math.sqrt((i - WaterSource.x)^2 + (j - WaterSource.y)^2)
			if dist <= radius then
				if i >= 0 and i <= GRID_WIDTH + 3 and j >= 0 and j <= GRID_HEIGHT + 3 then
					totalCells = totalCells + 1
				end
			end
		end
	end

	-- Now, add water to those cells
	for i = WaterSource.x - radius, WaterSource.x + radius do
		for j = WaterSource.y - radius, WaterSource.y + radius do
			local dist = math.sqrt((i - WaterSource.x)^2 + (j - WaterSource.y)^2)
			if dist <= radius then
				if i >= 0 and i <= GRID_WIDTH + 3 and j >= 0 and j <= GRID_HEIGHT + 3 then
					-- Adjust rate dynamically to simulate varying inflow (e.g., simulate rain)
					local dynamicRate = WaterSource.rate * (1 + 0.5 * math.sin(os.clock()))
					-- Add water to the cell
					H_Water[i][j] = H_Water[i][j] + (dynamicRate * deltaTime) / totalCells
					-- Update water surface elevation
					E[i][j] = H_Terrain[i][j] + H_Water[i][j]
				end
			end
		end
	end
end

-- Superbee Flux Limiter Function
local function Superbee(r)
	return math.max(0, math.max(math.min(2*r, 1), math.min(r, 2)))
end

-- Compute terrain slopes
function ComputeTerrainSlopes()
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			SlopeX[i][j] = (H_Terrain[i+1][j] - H_Terrain[i-1][j]) / (2 * DX)
			SlopeY[i][j] = (H_Terrain[i][j+1] - H_Terrain[i][j-1]) / (2 * DY)
		end
	end
end

-- Compute fluxes with flux limiting and viscosity
function ComputeFluxes(deltaTime)
	-- Compute terrain slopes first
	ComputeTerrainSlopes()

	-- Compute fluxes in x-direction
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 3 do
			local hL = H_Water[i-1][j]
			local hR = H_Water[i][j]
			local deltaHL = hR - hL
			local deltaHR = H_Water[i+1][j] - hR

			local r = deltaHL / (deltaHR + EPSILON)
			local phi = Superbee(r)

			local hAvg = math.min((hL + hR) / 2, MAX_WATER_HEIGHT)
			local deltaE = E[i][j] - E[i-1][j] + SlopeX[i][j] * DX

			-- Limit the slope
			local maxDeltaE = MAX_SLOPE * DX
			deltaE = math.clamp(deltaE, -maxDeltaE, maxDeltaE)

			local flux = -GRAVITY * deltaTime * hAvg * (deltaE / DX)
			flux = flux * phi  -- Apply flux limiter
			flux = flux * (1 - DAMPING_COEFFICIENT)  -- Apply damping

			-- Viscosity term
			local viscosityTerm = VISCOSITY_COEFFICIENT * (H_Water[i+1][j] - 2 * H_Water[i][j] + H_Water[i-1][j]) / (DX^2)
			flux = flux - viscosityTerm * deltaTime

			Qx[i][j] = flux
		end
	end

	-- Compute fluxes in y-direction
	for i = 1, GRID_WIDTH + 3 do
		for j = 1, GRID_HEIGHT + 2 do
			local hB = H_Water[i][j-1]
			local hT = H_Water[i][j]
			local deltaHB = hT - hB
			local deltaHT = H_Water[i][j+1] - hT

			local r = deltaHB / (deltaHT + EPSILON)
			local phi = Superbee(r)

			local hAvg = math.min((hB + hT) / 2, MAX_WATER_HEIGHT)
			local deltaE = E[i][j] - E[i][j-1] + SlopeY[i][j] * DY

			-- Limit the slope
			local maxDeltaE = MAX_SLOPE * DY
			deltaE = math.clamp(deltaE, -maxDeltaE, maxDeltaE)

			local flux = -GRAVITY * deltaTime * hAvg * (deltaE / DY)
			flux = flux * phi  -- Apply flux limiter
			flux = flux * (1 - DAMPING_COEFFICIENT)  -- Apply damping

			-- Viscosity term
			local viscosityTerm = VISCOSITY_COEFFICIENT * (H_Water[i][j+1] - 2 * H_Water[i][j] + H_Water[i][j-1]) / (DY^2)
			flux = flux - viscosityTerm * deltaTime

			Qy[i][j] = flux
		end
	end
end

-- Advect water heights using velocities
function Advect(deltaTime)
	local newH_Water = {}
	for i = 0, GRID_WIDTH + 3 do
		newH_Water[i] = {}
		for j = 0, GRID_HEIGHT + 3 do
			newH_Water[i][j] = H_Water[i][j]
		end
	end

	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			local u = (Qx[i][j] + Qx[i-1][j]) / (2 * (H_Water[i][j] + EPSILON))
			local v = (Qy[i][j] + Qy[i][j-1]) / (2 * (H_Water[i][j] + EPSILON))

			local dx = -u * deltaTime / DX
			local dy = -v * deltaTime / DY

			local i_upwind = i + dx
			local j_upwind = j + dy

			-- Bilinear interpolation (optional for better accuracy)
			local i0 = math.floor(i_upwind)
			local j0 = math.floor(j_upwind)
			local di = i_upwind - i0
			local dj = j_upwind - j0

			if i0 >= 0 and i0 <= GRID_WIDTH + 3 and j0 >= 0 and j0 <= GRID_HEIGHT + 3 then
				local h00 = H_Water[i0][j0]
				local h10 = H_Water[i0 + 1][j0]
				local h01 = H_Water[i0][j0 + 1]
				local h11 = H_Water[i0 + 1][j0 + 1]

				local h_interp = (1 - di) * (1 - dj) * h00 + di * (1 - dj) * h10 + (1 - di) * dj * h01 + di * dj * h11
				newH_Water[i][j] = h_interp
			end
		end
	end

	-- Update water heights
	H_Water = newH_Water
end

-- Update water heights and apply numerical diffusion
function UpdateHeights()
	local newH_Water = {}
	for i = 0, GRID_WIDTH + 3 do
		newH_Water[i] = {}
		for j = 0, GRID_HEIGHT + 3 do
			newH_Water[i][j] = H_Water[i][j]
		end
	end

	-- Update water heights based on fluxes
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			local fluxInX = Qx[i-1][j]
			local fluxOutX = Qx[i][j]
			local fluxInY = Qy[i][j-1]
			local fluxOutY = Qy[i][j]

			local netFlux = (fluxInX - fluxOutX) / DX + (fluxInY - fluxOutY) / DY
			newH_Water[i][j] = H_Water[i][j] + netFlux

			-- Prevent negative water height and handle wet/dry boundaries
			if newH_Water[i][j] < MIN_WATER_DEPTH then
				newH_Water[i][j] = 0
				Qx[i][j] = 0
				Qy[i][j] = 0
			elseif newH_Water[i][j] > MAX_WATER_HEIGHT then
				-- Limit maximum water height
				newH_Water[i][j] = MAX_WATER_HEIGHT
			end
		end
	end

	-- Apply advanced diffusion handling
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			local u = (Qx[i][j] + Qx[i-1][j]) / (2 * (H_Water[i][j] + EPSILON))
			local v = (Qy[i][j] + Qy[i][j-1]) / (2 * (H_Water[i][j] + EPSILON))
			local speed = math.sqrt(u^2 + v^2)

			local localDiffusion = DIFFUSION_COEFFICIENT * (1 + speed)
			local laplacian = H_Water[i-1][j] + H_Water[i+1][j] + H_Water[i][j-1] + H_Water[i][j+1] - 4 * H_Water[i][j]
			newH_Water[i][j] = newH_Water[i][j] + localDiffusion * laplacian

			-- Ensure non-negative water height
			if newH_Water[i][j] < MIN_WATER_DEPTH then
				newH_Water[i][j] = 0
			end
		end
	end

	-- Update water heights and surface elevations
	for i = 0, GRID_WIDTH + 3 do
		for j = 0, GRID_HEIGHT + 3 do
			H_Water[i][j] = newH_Water[i][j]
			E[i][j] = H_Terrain[i][j] + H_Water[i][j]
		end
	end
end

-- Apply smoothing to water heights
function SmoothWaterHeights()
	local smoothedH_Water = {}
	for i = 0, GRID_WIDTH + 3 do
		smoothedH_Water[i] = {}
		for j = 0, GRID_HEIGHT + 3 do
			smoothedH_Water[i][j] = H_Water[i][j]
		end
	end

	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			smoothedH_Water[i][j] = (H_Water[i][j] + H_Water[i-1][j] + H_Water[i+1][j] + H_Water[i][j-1] + H_Water[i][j+1]) / 5
		end
	end

	H_Water = smoothedH_Water
end

-- Enforce boundary conditions
function EnforceBoundaryConditions()
	-- Absorbing boundary conditions
	for i = 0, GRID_WIDTH + 3 do
		H_Water[i][0] = H_Water[i][1] * 0.9
		H_Water[i][GRID_HEIGHT + 3] = H_Water[i][GRID_HEIGHT + 2] * 0.9
		Qy[i][0] = Qy[i][1] * 0.9
		Qy[i][GRID_HEIGHT + 3] = Qy[i][GRID_HEIGHT + 2] * 0.9
	end
	for j = 0, GRID_HEIGHT + 3 do
		H_Water[0][j] = H_Water[1][j] * 0.9
		H_Water[GRID_WIDTH + 3][j] = H_Water[GRID_WIDTH + 2][j] * 0.9
		Qx[0][j] = Qx[1][j] * 0.9
		Qx[GRID_WIDTH + 3][j] = Qx[GRID_WIDTH + 2][j] * 0.9
	end
end

-- Check and enforce water conservation
function ConserveWater(previousTotalWater)
	local currentTotalWater = 0
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			currentTotalWater = currentTotalWater + H_Water[i][j] * DX * DY
		end
	end
	-- Account for water added this step
	-- Water conservation correction can be applied if necessary
end

-- Rendering
function InitializeWaterParts()
	local WaterModel = Instance.new("Model")
	WaterModel.Name = "Water"
	WaterModel.Parent = workspace

	for i = 0, GRID_WIDTH + 3 do
		WaterParts[i] = {}
		for j = 0, GRID_HEIGHT + 3 do
			local part = Instance.new("Part")
			part.Size = Vector3.new(DX, 1, DY)
			part.Anchored = true
			part.CanCollide = false
			part.Material = Enum.Material.Water
			part.Color = Color3.fromRGB(0, 162, 255)
			part.Transparency = 0.5
			part.Parent = WaterModel
			WaterParts[i][j] = part
		end
	end
end

function RenderSimulation()
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			local waterDepth = H_Water[i][j]
			local terrainHeight = H_Terrain[i][j]
			local part = WaterParts[i][j]

			if waterDepth > MIN_WATER_DEPTH then
				-- Position the water part
				part.Position = Vector3.new(
					(i - 2) * DX + OFFSET_X,
					terrainHeight + waterDepth / 2,
					(j - 2) * DY + OFFSET_Z
				)
				-- Adjust the size (height)
				part.Size = Vector3.new(DX, waterDepth, DY)
				-- Adjust color based on depth
				local colorIntensity = math.clamp(waterDepth / MAX_WATER_HEIGHT, 0, 1)
				part.Color = Color3.fromRGB(0, math.floor(162 * colorIntensity), 255)
				part.Transparency = 0.5
			else
				-- Hide the part if no water
				part.Size = Vector3.new(0, 0, 0)
				part.Transparency = 1
			end
		end
	end
end

-- Initialize the simulation
GenerateTerrainHeightmap()
InitializeWaterParts()

-- Initial total water volume
local totalWaterVolume = 0
for i = 1, GRID_WIDTH + 2 do
	for j = 1, GRID_HEIGHT + 2 do
		totalWaterVolume = totalWaterVolume + H_Water[i][j] * DX * DY
	end
end

-- Main simulation loop
while true do
	-- Compute maximum water depth for time step calculation
	local maxWaterDepth = MIN_WATER_DEPTH
	for i = 1, GRID_WIDTH + 2 do
		for j = 1, GRID_HEIGHT + 2 do
			if H_Water[i][j] > maxWaterDepth then
				maxWaterDepth = H_Water[i][j]
			end
		end
	end

	-- Compute time step based on CFL condition
	local deltaTime = CFL_NUMBER * math.min(DX, DY) / (math.sqrt(GRAVITY * maxWaterDepth) + EPSILON)
	deltaTime = math.min(deltaTime, 0.01)  -- Cap deltaTime to a safe maximum

	-- Simulation steps
	EnforceBoundaryConditions()
	ComputeFluxes(deltaTime)
	Advect(deltaTime)
	UpdateHeights()
	SmoothWaterHeights()
	AddWater(deltaTime)
	ConserveWater(totalWaterVolume)
	RenderSimulation()

	-- Wait for the next time step
	wait(deltaTime)
end

Currently it goes crazy and dissapears when it starts and somethimes freezes and throws the error ServerScriptService.4:250: attempt to perform arithmetic (mul) on number and nil  -  Server - 4:250
				local h_interp = (1 - di) * (1 - dj) * h00 + di * (1 - dj) * h10 + (1 - di) * dj * h01 + di * dj * h11