SandDunesConfig.lua

--------------------------------------------------------------------------------
-- Initial Parameters
--------------------------------------------------------------------------------
local myProblem = {}

myProblem.Init = function(self, o)

		-- Numerical parameters of the discretization
		
	self.problem = o.problem
	self.simCase = o.simCase
	self.simCaseBnd = o.simCaseBnd
	
	self.dim = o.dim
	self.dir_name = o.dir_name
	self.file_name = o.file_name
	self.folder_name = o.folder_name
	self.elem_type = o.elem_type
	self.numRefs = o.numRefs
	self.numPreRefs = o.numPreRefs
	self.algebraBlockSize = o.algebraBlockSize
	
	self.timeMethod  = o.timeMethod
	self.DT = o.DT
	self.DTmax = o.DTmax
	self.DTmin = o.DTmin
	self.DTLimex = o.DTLimex
	self.numTimeSteps = o.numTimeSteps
	self.startTime = o.startTime
	self.endTime = o.endTime
	self.modifyDT  = o.modifyDT
	
	
		--Output Data
	self.boolData = o.boolData
	self.data_name = o.data_name
	self.outputFactor = o.outputFactor
	self.writeIntegral = o.writeIntegral
	self.boolLoadCheckPoint = o.boolLoadCheckPoint
	self.boolSaveCheckPoint = o.boolSaveCheckPoint

	self.incr_factor = o.incr_factor
	self.red_factor_fail = o.red_factor_fail
	self.red_factor_success = o.red_factor_success
	self.optimal_newton_steps = o.optimal_newton_steps
	self.maxConvRate = o.maxConvRate
	self.minConvRate = o.minConvRate
	self.NewtonDebug = o.NewtonDebug
	self.NewtonSteadyDebug = o.NewtonSteadyDebug
	self.NewtonUpdaterSteady = o.NewtonUpdaterSteady
	self.NewtonUpdater = o.NewtonUpdater
	self.StepDebug = o.StepDebug
	self.debug_dir = o.debug_dir
	
		
	self.tol = o.tol
	self.nstages = o.nstages
	self.limex_partial_mask = o.limex_partial_mask
	self.limex_debug_level = o.limex_debug_level
	self.VelErrorNorm = o.VelErrorNorm
	self.PressErrorNorm = o.PressErrorNorm
	self.VolErrorNorm = o.VolErrorNorm
	self.alphaVel  = o.alphaVel
	self.alphaPress = o.alphaPress
	self.alphaVol = o.alphaVol

	self.max_newton_steps_steady_state = o.max_newton_steps_steady_state
	self.max_newton_steps_transient = o.max_newton_steps_transient
	
	self.SteadyAbsDefect = o.SteadyAbsDefect
	self.SteadyRedDefect = o.SteadyRedDefect
	
	self.AbsDefect = o.AbsDefect
	self.RedDefect = o.RedDefect

	self.baseSolverLim = o.baseSolverLim
	self.baseSolverImp = o.baseSolverImp
	
	self.LinAbsDefectImp = o.LinAbsDefectImp
	self.LinRedDefectImp = o.LinRedDefectImp
	self.LinAbsDefectLim = o.LinAbsDefectLim
	self.LinRedDefectLim = o.LinRedDefectLim
	self.max_linear_steps_Lim = o.max_linear_steps_Lim
	self.max_linear_steps_Imp = o.max_linear_steps_Imp
	self.rap = o.rap
	self.rim = o.rim
		
	self.damping_gmg = o.damping_gmg
	self.damping_ilu = o.damping_ilu
	self.value_beta = o.value_beta
	
	self.linear_solver_Lim = o.linear_solver_Lim
	self.precondLim = o.precondLim
	self.smootherLim = o.smootherLim
	self.orderingLim = o.orderingLim
	self.pre_smoothLim = o.pre_smoothLim
	self.post_smoothLim = o.post_smoothLim
	self.gmgBaseLevelLim = o.gmgBaseLevelLim

	self.linear_solver_Imp = o.linear_solver_Imp
	self.precondImp = o.precondImp
	self.smootherImp = o.smootherImp
	self.orderingImp = o.orderingImp
	self.pre_smoothImp = o.pre_smoothImp
	self.post_smoothImp = o.post_smoothImp
	self.gmgBaseLevelImp = o.gmgBaseLevelImp

	
	self.eps_ilut = o.eps_ilut

	self.lambdamaxSteps = o.lambdamaxSteps
	self.lambdaStart = o.lambdaStart
	
		-- Physical phenomenon of simulation
	self.doSteadyState = o.doSteadyState
	self.steadyStateOnly = o.steadyStateOnly
	self.boolSource = o.boolSource
	self.consistentRho_in_source = o.consistentRho_in_source
	self.boolRelativeVel = o.boolRelativeVel
	self.boolGradientPsSource = o.boolGradientPsSource
	self.boolViscPs = o.boolViscPs
	self.boolAveDiff = o.boolAveDiff
	self.boolSaltationFlux = o.boolSaltationFlux
	self.boolSlipDiff = o.boolSlipDiff
	self.boolSlipVel = o.boolSlipVel
	self.boolpress_jump = o.boolpress_jump
	self.boolAveNormal = o.boolAveNormal
	self.boolFixVel = o.boolFixVel
	self.boolFixVol = o.boolFixVol
	self.boolMassTerm = o.boolMassTerm
	self.boolDensityMean = o.boolDensityMean
	
	self.inflow = o.inflow
	self.SlipVelValue = o.SlipVelValue
	self.H_0 = o.H_0
	self.ReferencePressure = o.ReferencePressure
	self.bStokes = o.bStokes or false
	self.bNoLaplace = o.bNoLaplace or false
	self.exactJac = o.exactJac or false
	self.bPecletBlend = o.bPecletBlend or false
	self.upwind_m = o.upwind_m or "full"
	self.upwind_t = o.upwind_t or "full"
	self.upwind_r = o.upwind_r or "full"
	self.bPac = o.bPac
	self.diffLength = o.diffLength
	self.stab = o.stab
	self.div_correction = o.div_correction
	self.boolIPVelocity = o.boolIPVelocity
	self.turbViscMethod = o.turbViscMethod
	self.modellconstant = o.modellconstant
	self.beta1 = o.beta1
	self.kIn = o.kIn
	self.omegaIn = o.omegaIn
	self.dWall = o.dWall
	self.omegaWall = o.omegaWall

	
	--Material Properties
	self.nu_a = o.nu_a
	self.rho_a = o.rho_a
	self.rho_s = o.rho_s
	self.dp = o.dp
	self.nu_s = o.nu_s
	self.c_init = o.c_init

	self.alpha_max = o.alpha_max
	self.alpha_min = o.alpha_min
	self.packing_factor = o.packing_factor
	self.grad_limit = o.grad_limit
	self.slope_limit = o.slope_limit
	self.granular_model = o.granular_model
	self.density_model = o.density_model
	self.interface_value = o.interface_value
	self.drag_mod = o.drag_mod
	self.riemman = o.riemman


	self.FR = o.FR
	self.B_phi = o.B_phi
	self.deltaGamma = o.deltaGamma
	self.Visc_limit = o.Visc_limit

	self.deltaPs = o.deltaPs
	self.deltaI = o.deltaI
	self.FricMu_1 = o.FricMu_1
	self.FricMu_2 = o.FricMu_2
	self.I_0 = o.I_0
	self.gravity = o.gravity
	self.roughness_length = o.roughness_length

	
	
	self.boolSolverDesc = false
	self.NewtonSolverDescSteady = {}
	self.NewtonSolverDesc = {}
	
	
	self.gridName = o.gridName
  
  
end
--------------------------------------------------------------------------------
-- DEFAULT Parameters
--------------------------------------------------------------------------------
myProblem.GetCaseDefaults = function(self, problem)
	if problem == "flow" then
		return {
			Name = "Dune",
			numRefs = 6,
			numPreRefs = 3,
			alphaVel = 1.0,
			alphaPress = 1e-11,
			boolSaltationFlux = false,
			boolAveDiff = false,
			boolSlipDiff = false,
			boolSlipVel = false,
			boolFixVel = false,
			rap = true,
			rim = false,
			value_beta = 0.0,
			
			linear_solver_Imp = "bicgStab",
			smootherImp = "ilu",
			orderingImp = "none",
			pre_smoothImp = 1,
			post_smoothImp = 0,
			baseSolverImp = "SuperLu",
			LinAbsDefectImp = 1e-08,
			LinRedDefectImp = 1e-3,
			gmgBaseLevelImp = 5,
			
			linear_solver_Lim = "gmres",
			smootherLim = "ilu",
			orderingLim = "none",
			pre_smoothLim = 1,
			post_smoothLim = 0,
			baseSolverLim = "SuperLu",
			LinAbsDefectLim = 1e-8,
			LinRedDefectLim = 1e-7,
			gmgBaseLevelLim = 5,
			
			DTLimex=1000,
		}
	elseif problem == "avalanche" then
		return {
			Name = "Solution",
			numRefs = 4,
			numPreRefs = 3,
			alphaVel = 1e-15,
			alphaPress = 0.5e-15,
			boolSaltationFlux = false,
			boolAveDiff = false,
			boolSlipDiff = false,
			boolSlipVel = true,
			boolFixVel = true,
			rap = true,
			rim = false,
			value_beta = 0.0,
			
			linear_solver_Imp = "bicgStab",
			smootherImp = "ilu",
			orderingImp = "none",
			pre_smoothImp = 1,
			post_smoothImp = 2,
			baseSolverImp = "SuperLu",
			LinAbsDefectImp = 1e-08,
			LinRedDefectImp = 1e-5,
			gmgBaseLevelImp = 5,
			
			linear_solver_Lim = "bicgStab",
			smootherLim = "ilu",
			orderingLim = "none",
			pre_smoothLim = 1,
			post_smoothLim = 2,
			baseSolverLim = "SuperLu",
			LinAbsDefectLim = 1e-12,
			LinRedDefectLim = 1e-12,
			gmgBaseLevelLim = 5,
			
			
			DTLimex=1e-03,
		}
	end

	error("Unknown problem: " .. tostring(problem))
end
--------------------------------------------------------------------------------
-- GEOMETRY
--------------------------------------------------------------------------------
myProblem.GetGeometry = function(self, params)
    if params.dim ~= 2 and params.dim ~= 3 then error("Only dimensions 2 and 3 are supported") end
    if params.elem_type ~= "tri" and params.elem_type ~= "quad" then error("Unknown element type: " .. tostring(params.elem_type)) end

    local geometry = {}

    if params.problem == "flow" then
        geometry.gridName = "Dune" .. params.dim .. "D_" .. params.elem_type .. "_double.ugx"
        geometry.innerSubsets = {"Inner", "Inner2"}
        geometry.allSubsets = "Inner,Inner2,Left,Right,Top,Bottom"
        geometry.boundaries = "Left,Right,Top,Bottom"

        if params.dim == 3 then
            geometry.allSubsets = geometry.allSubsets .. ",Back,Front"
        end
    elseif params.problem == "avalanche" then
        geometry.gridName = "Avalanche" .. params.dim .. "D_" .. params.elem_type .. ".ugx"
        geometry.innerSubsets = {"Inner"}

        if params.dim == 2 then
            geometry.allSubsets = "Inner,Left,Right,Top,Bottom"
            geometry.boundaries = "Left,Right,Top,Bottom"
        else
            geometry.allSubsets = "Inner,Left1,Left2,Right,Top,Bottom,Front1,Front2,Back1,Back2,Back3"
            geometry.boundaries = "Left1,Left2,Right,Top,Bottom,Front1,Front2,Back1,Back2,Back3"
        end
    else
        error("Unknown problem: " .. tostring(params.problem))
    end
	self.geometry = geometry
    return geometry
end
--------------------------------------------------------------------------------
-- File Names
--------------------------------------------------------------------------------
myProblem.FileNames = function (self,rank,SpaceSize)


	local riemman_name = nil
	if self.riemman == 0 then
		riemman_name = "Godunov"
	elseif self.riemman == 1 then
		riemman_name = "Rusanov"
	elseif self.riemman == 2 then
		riemman_name = "Roe"
	elseif self.riemman == 3 then
		riemman_name = "Engquist-Osher"
	else
		print ("Numerical Flux Scheme for NonLinear Scalar conservation Law not defined"); exit();
	end

	local folder_name = nil
	
	if self.dir_name == "" then
		folder_name = self.folder_name
	else
		folder_name = self.dir_name .. "/" .. self.folder_name
	end

	---------------------
	-- FolderName
	---------------------

	folder_name = folder_name .. "_" .. self.dim.. "D"
	folder_name = folder_name .. "-" .. self.timeMethod
	if self.turbViscMethod ~= "no" and not(self.bStokes) then
		folder_name = folder_name .. "-" .. self.turbViscMethod
	end
	if self.bStokes then
		folder_name = folder_name .. "-Stokes"
	end

	folder_name = folder_name .. "-" ..riemman_name
	
	if self.boolSaltationFlux then
		folder_name =folder_name .. "-SaltFlux"
	else
		folder_name =folder_name .. "-NoSaltFlux"
	end

	if self.boolRelativeVel then
		folder_name =folder_name .. "-RelVel"
	else
		folder_name =folder_name .. "-NoRelVel"
	end


	if  self.boolGradientPsSource  or self.boolSource  then

		if self.boolSource then
			if self.consistentRho_in_source then
				folder_name =folder_name .. "-Consistent"
			end
			folder_name = folder_name .. "-MG_Force"
		end
		
		if self.boolGradientPsSource then
			folder_name = folder_name .. "-DPs"
		else
			folder_name = folder_name .. "-DPh"
		end
		
	else
		folder_name = folder_name .. "-NoForce"
	end

	if self.boolAveDiff then
		folder_name = folder_name .. "-AveDiff"
	else
		folder_name = folder_name .. "-NoAveDiff"
	end

	if (self.boolSlipDiff and self.boolSlipVel) then print ("SlipVel and SlipDiff activated at the same time."); exit(); end
	if self.boolSlipDiff then
		folder_name = folder_name .. "-SlipDiff"
	else if self.boolSlipVel then
			folder_name = folder_name .. "-SlipVel"
		else
			folder_name = folder_name .. "-NoSlip"
		end
	end
	

	if (rank==0) then
		if not DirectoryExists (folder_name) then
			CreateDirectory (folder_name)
		end
		print(folder_name)
	end
	
	---------------------
	-- VTK folder
	---------------------
	local file_name = self.file_name
	
	vtk_file_name = file_name .. "-" .. self.elem_type
	vtk_file_name = vtk_file_name .. "-lev" .. self.numRefs
	
	local fixedSpaceSize = string.format("%04d", SpaceSize)
	vtk_file_name = vtk_file_name .. "_Parallel_" .. fixedSpaceSize

	folder_vtk = folder_name .. "/" .. vtk_file_name
	
	if (rank== 0) then
		
		if not DirectoryExists (folder_vtk) then
			CreateDirectory (folder_vtk)
		else
			if(self.NewtonDebug or self.NewtonSteadyDebug) then
				local cmd = "rm -rf " .. folder_vtk
				os.execute(cmd)
				print(" Directory reseted")
				CreateDirectory (folder_vtk)
			end
		end
	end

	

	vtk_file_name = folder_vtk .. "/" .. self.file_name -- VTK output file name base
	
	
	self.debug_dir = folder_vtk
	self.riemman_name = riemman_name

	print("Files Setting DONE")
	print("File Adress: ".. vtk_file_name)
	
	-- vtk_file_name: Name of the vtk files
	-- folder_vtk: Directory where vtk files are saved
	-- folder_name: Directory where the different simulations are saved, e.g. folder2dLimex, folder3dEuler
	
	return vtk_file_name,folder_vtk,folder_name
end

--------------------------------------------------------------------------------
-- LOG Files
--------------------------------------------------------------------------------
myProblem.LogFiles = function (self,rank_t,Name)
	local counter = 0
	local filename = Name .. "_" .. counter

	while io.open(filename, "r") do
		counter = counter + 1
		filename = Name .. "_" .. counter
	end

	GetLogAssistant():enable_file_output(true, filename)

	if rank_t > 0 then
		GetLogAssistant():enable_terminal_output(false)
	end
end





--------------------------------------------------------------------------------
-- Printing Simulation Settings
--------------------------------------------------------------------------------
myProblem.PrintingSettings = function (self)
	-----------------------------------------------------------------------
	print (" Sand Dune Dynamics     " .. os.date("%A, %B %d, %Y at %I:%M %p"))
	print (" Geometry: " .. self.gridName ..", dim = " .. self.dim)
	print("Algebra Block Size = " .. tostring(self.algebraBlockSize))
	print (" Physical parameter:")
	print ("	Table case		= " .. self.simCase + 1)
	print ("	inflow			= " .. self.inflow)
	print ("	SlipVel			= " .. self.SlipVelValue)
	print ("	Stokes			= " .. tostring (self.bStokes))
	print ("	Steady state    	= " .. tostring (self.doSteadyState))
	print ("	BodyForce       	= " .. tostring (self.boolSource))
	print ("	Consisten Gravity	= " .. tostring (self.consistentRho_in_source))
	print ("	Saltation FLux       	= " .. tostring (self.boolSaltationFlux))
	print ("	Relative Vel    	= " .. tostring (self.boolRelativeVel))
	print ("	Ps gradient     	= " .. tostring (self.boolGradientPsSource))
	print ("	Ps in visc      	= " .. tostring (self.boolViscPs))
	print ("	Diffusion       	= " .. tostring (self.boolAveDiff))
	print ("	SlipDiff         	= " .. tostring (self.boolSlipDiff))
	print ("	SlipVel         	= " .. tostring (self.boolSlipVel))
	print ("	MassMean         	= " .. tostring (self.boolDensityMean))
	-----------------------------------------------------------------------
	print (" Numerical parameter:")
	print ("	elem_type		= " .. self.elem_type)
	print ("	numRefs			= " .. self.numRefs)
	print ("	numPreRefs		= " .. self.numPreRefs)
	print ("	timeMethod		= " .. self.timeMethod)
	print ("	DT   			= " .. self.DT)
	print ("	numTimeSteps		= " .. self.numTimeSteps)
	print ("	noLaplace		= " .. tostring (self.bNoLaplace))
	print ("	exactJac		= " .. tostring (self.exactJac))
	print ("	PecletBlend		= " .. tostring (self.bPecletBlend))
	print ("	upwind_m		= " .. self.upwind_m)
	print ("	upwind_t		= " .. self.upwind_t)
	print ("	Num Flux Scheme		= " .. self.riemman_name)
	print ("	pac			= " .. tostring (self.bPac))
	print ("	stab			= " .. self.stab)
	print ("	difflength		= " .. self.diffLength)
	print ("	Turbulence		= " .. self.turbViscMethod)
	-----------------------------------------------------------------------
	print (" Linear Solver parameters:")
	print ("	damping_gmg		= " .. self.damping_gmg)
	print ("	damping_ilu		= " .. self.damping_ilu)
	print ("	RAP			= " .. tostring (self.rap))
	print ("	RIM			= " .. tostring (self.rim))
	print ("	beta			= " .. self.value_beta)
	print ("	AbsDefImp		= " .. self.LinAbsDefectImp)
	print ("	RedDefImp		= " .. self.LinRedDefectImp)
	print ("	AbsDefLim 		= " .. self.LinAbsDefectLim)
	print ("	RedDefLim		= " .. self.LinRedDefectLim)
	print ("	MaxStepsLim		= " .. tostring (self.max_linear_steps_Lim))
	print ("	MaxStepsImp		= " .. tostring (self.max_linear_steps_Imp))
	
	print ("	PreconditionerLim	= " .. tostring (self.precondLim))
	print ("	Smoother Lim		= " .. tostring (self.smootherLim))
	print ("	Ordering Lim		= " .. tostring (self.orderingLim))
	print ("	Pre Smooth Lim		= " .. tostring (self.pre_smoothLim))
	print ("	Post Smooth Lim		= " .. tostring (self.post_smoothLim))
	
	print ("	PreconditionerImp	= " .. tostring (self.precondImp))
	print ("	Smoother Imp		= " .. tostring (self.smootherImp))
	print ("	Ordering Imp		= " .. tostring (self.orderingImp))
	print ("	Pre Smooth Imp		= " .. tostring (self.pre_smoothImp))
	print ("	Post Smooth Imp		= " .. tostring (self.post_smoothImp))
	print ("	Eps ILUT		= " .. self.eps_ilut)
	
	-----------------------------------------------------------------------

	if self.timeMethod == "limex" then
		print (" Limex Numerical parameter:")
		print ("	Vel norm		= " .. self.VelErrorNorm)
		print ("	Press norm		= " .. self.PressErrorNorm)
		print ("	Vol  norm		= " .. self.VolErrorNorm)
		print ("	Vel norm scale		= " .. self.alphaVel)
		print ("	Press norm scale	= " .. self.alphaPress)
		print ("	Vol norm scale		= " .. self.alphaVol)
	end

end


--------------------------------------------------------------------------------
-- Approximation Space
--------------------------------------------------------------------------------
myProblem.ApproximationSpace = function (self)
	local allSubsets = self.geometry.allSubsets
	print("Approximation Space Setting")
	local fct_cmp_tbl = nil
	local nv_cmp_tbl = nil
	local vel_cmp_tbl = nil
	if self.turbViscMethod == "komegaSST" then
		if self.dim == 3 then
			fct_cmp_tbl = {"u", "v", "w", "p", "c", "k", "omega"}
			nv_cmp_tbl = {"u", "v", "w", "p", "c"}
			vel_cmp_tbl = {"u", "v", "w"}
		else
			fct_cmp_tbl = {"u", "v", "p", "c", "k", "omega"}
			nv_cmp_tbl = {"u", "v", "p", "c"}
			vel_cmp_tbl = {"u", "v"}
		end
	else
		if self.dim == 3 then
			fct_cmp_tbl = {"u", "v", "w", "p", "c"}
			nv_cmp_tbl = fct_cmp_tbl
			vel_cmp_tbl = {"u", "v", "w"}
		else
			fct_cmp_tbl = {"u", "v", "p", "c"}
			nv_cmp_tbl = fct_cmp_tbl
			vel_cmp_tbl = {"u", "v"}
		end
	end

	-- Create the domain, load the grid and refine it
	local dom = util.CreateDomain (self.gridName, self.numPreRefs)
	balancer.RefineAndRebalanceDomain (dom, self.numRefs - self.numPreRefs)
	--util.refinement.CreateRegularHierarchy(dom, self.numRefs, true)

	print ("Domain-info:")
	print (dom:domain_info():to_string())

	-- Create the vertex-centered approximation space
	local approxSpace = ApproximationSpace (dom)

	approxSpace:add_fct("u", "Lagrange",1,allSubsets)
	approxSpace:add_fct("v", "Lagrange",1,allSubsets)
	if self.dim == 3 then
		approxSpace:add_fct("w", "Lagrange",1,allSubsets)
	end
	approxSpace:add_fct("p", "Lagrange",1,allSubsets)
	approxSpace:add_fct("c", "Lagrange",1,allSubsets)
	
	if self.turbViscMethod == "komegaSST" then
		approxSpace:add_fct("k", "Lagrange",1,allSubsets)
		approxSpace:add_fct("omega", "Lagrange",1,allSubsets)
	end

	approxSpace:init_levels()
	approxSpace:init_top_surface()
	approxSpace:print_statistic()

	OrderLex (approxSpace, "x")
	--OrderCuthillMcKee(approxSpace,true)

	util.solver.defaults.approxSpace = approxSpace


	-- grid function for the solution
	local u = GridFunction (approxSpace)
	u:set(0)

	self.fct_cmp_tbl = fct_cmp_tbl
	self.nv_cmp_tbl = nv_cmp_tbl
	self.vel_cmp_tbl = vel_cmp_tbl
	self.u = u
	self.approxSpace = approxSpace
	
	print("Approximation Space DONE")
	
	
	return approxSpace,u
end

--------------------------------------------------------------------------------
-- Interface Parameters List
--------------------------------------------------------------------------------
myProblem.InterfaceParameters = function (self)
	
	local InterfaceValues = Interface()

	InterfaceValues:set_particle_diameter(self.dp)
	InterfaceValues:set_particle_density(self.rho_s)
	InterfaceValues:set_air_density(self.rho_a)
	InterfaceValues:set_fluid_Visc(self.nu_a*self.rho_a)
	InterfaceValues:set_particle_kinVisc(self.nu_s)
	InterfaceValues:set_alpha_max(self.alpha_max)
	InterfaceValues:set_alpha_min(self.alpha_min)
	InterfaceValues:set_packing_factor(self.packing_factor)
	InterfaceValues:set_FrictionMu_1(self.FricMu_1)
	InterfaceValues:set_FrictionMu_2(self.FricMu_2)
	InterfaceValues:set_I_0(self.I_0)
	InterfaceValues:set_deltaI(self.deltaI)
	InterfaceValues:set_deltaPs(self.deltaPs)
	InterfaceValues:set_FR(self.FR)
	InterfaceValues:set_B_phi(self.B_phi)
	InterfaceValues:set_deltaGamma(self.deltaGamma)
	InterfaceValues:set_limit(self.Visc_limit)
	InterfaceValues:set_bool_particle_pressure_force(self.boolGradientPsSource)
	InterfaceValues:set_bool_consistent_gravity(false)
	InterfaceValues:set_reference_pressure(self.ReferencePressure)
	InterfaceValues:set_time_step_factor(self.DT)
	InterfaceValues:set_interface_volume_fraction(self.interface_value)
	InterfaceValues:set_drag_model(self.drag_mod)
	InterfaceValues:set_relative_vel_error(1)
	InterfaceValues:set_bool_initialized(true)
	
	self.InterfaceValues = InterfaceValues
	
	print("Interfce Parameters list: DONE")
	
	return InterfaceValues
end


--------------------------------------------------------------------------------
-- LUA Functions
--------------------------------------------------------------------------------
myProblem.RegisterCallbacks = function(self)
    if self.problem == "avalanche" then
        if self.simCaseBnd ~= 1 and self.simCaseBnd ~= 2 then error("Unknown avalanche boundary case") end
        if self.dim == 3 and self.simCaseBnd == 2 then error("The supplied Avalanche.lua has no initial concentration function for 3D case 2") end
    end

    local function inletSpeed(height)
		local hh=20.0
        if self.problem == "avalanche" then return 0.0 end
		if height <= 0.0 then return 0.0 end
		--return self.inflow
        return self.inflow * math.log(1 + height / self.roughness_length) / math.log(1 + hh / self.roughness_length)
		--return self.inflow* (2*hh - height) * (height ) / (hh * hh)
    end
	local function InflowVel2d(x, y, t) return inletSpeed(y), 0.0 end
	local function InflowVel3d(x, y, z, t) return inletSpeed(z), 0.0, 0.0 end

    local function initialC2d(x, y)
		local sigma = 5
        if self.problem == "flow" then
            local dx = math.abs((x - 20.0) / sigma)
            local height = self.H_0 * math.exp(-0.5*(dx^4)) - 0.01
            if y > height then return 0.0 end
            return self.c_init
        end

        if (x > 0.2 and x < 0.4) or x > 0.8 then return 1.5e-3 end
        return 0.0
    end

    local function initialC3d(x, y, z)
		local sigma = 5
        if self.problem == "avalanche" then return 0.0 end
        local dx = math.abs((x - 20.0) / sigma)^2
        local dy = math.abs(y / sigma)^2
        local height = self.H_0 * math.exp(-0.5*(dx + dy)^2) - 0.01
        if z > height then return 0.0 end
        return self.c_init
    end

    local function topFlux2d(x, y)
        if self.simCaseBnd == 1 then
            if x > 0.3 and x < 2.0 then return -1e-5 end
        else
            if (x > 0.2 and x < 0.4) or x > 0.8 then return -1e-5 end
        end
        return 0.0
    end

    local function topFlux3d(x, y, z)
        if x > -1.0 and x < 2.0 then return -1e-5 end
        return 0.0
    end
    
	local turbulenceIntensity = 0.05
	local turbulenceLength = 0.05
	local Cmu = 0.09
	
	local function InletK2d(x, y, t)
		local ux, uy = InflowVel2d(x, y, t)
		local U = math.sqrt(ux * ux + uy * uy)
		return 1.5 * (turbulenceIntensity * U)^2
	end
	local function InletK3d(x, y, z, t)
		local ux, uy, uz = InflowVel3d(x, y, z, t)
		local U = math.sqrt(ux * ux + uy * uy + uz * uz)
		return 1.5 * (turbulenceIntensity * U)^2
	end
	local function InletOmega2d(x, y, t)
		local k = InletK2d(x, y, t)
		return math.sqrt(k) / (Cmu^0.25 * turbulenceLength)
	end
	local function InletOmega3d(x, y, z, t)
		local k = InletK3d(x, y, z, t)
		return math.sqrt(k) / (Cmu^0.25 * turbulenceLength)
	end
	local function WallDistance2d(x, y, t)
		local lowerWallY = 0.0
		local dLower = y - lowerWallY
		return math.max( dLower, 0.0)
	end
	local function WallDistance3d(x, y, z, t)
		local lowerWallZ = 0.0
		local dLower = z - lowerWallZ
		return math.max( dLower, 0.0)
	end

    if self.dim == 2 then
        _G.SandInflow = InflowVel2d
        _G.SandInitialC = initialC2d
		_G.SandInitialK = InletK2d
		_G.SandInitialOmega = InletOmega2d
        _G.SandTopFlux = topFlux2d
		_G.InletK = InletK2d
		_G.InletOmega = InletOmega2d
		_G.WallDistance = WallDistance2d
    else
        _G.SandInflow = InflowVel3d
        _G.SandInitialC = initialC3d
		_G.SandInitialK = InletK3d
		_G.SandInitialOmega = InletOmega3d
        _G.SandTopFlux = topFlux3d
		_G.InletK = InletK3d
		_G.InletOmega = InletOmega3d
		_G.WallDistance = WallDistance3d
    end
end
--------------------------------------------------------------------------------
-- Variables
--------------------------------------------------------------------------------
myProblem.Clousures = function (self,u,Inner_total)
		
	local InterfaceValues = self.InterfaceValues
	-------------------------------------------------------------- VelocityGradMag

	--local gamma = ShearStressFV1(self.approxSpace,u)
	--self.gamma = gamma

	---------------------------------------------------------------------- Density
	local Density = nil
	if self.bStokes or self.boolSaltationFlux then
		Density = ConstUserNumber(self.rho_a)
	else
		Density = GranularDensityLinker();
		Density:set_model(self.density_model)
		Density:set_phase_parameters(self.InterfaceValues)
	end


	self.Density = Density

	---------------------------------------------------------------------- Viscosity

	local Inverse_RHO = InverseLinker();
	Inverse_RHO:divide(1.0,Density);

	local AirDensity = ConstUserNumber(self.rho_a)

	local Scale_RHO = ScaleLinker();
	Scale_RHO:set_import_1(Inverse_RHO)
	Scale_RHO:set_import_2(AirDensity)

	local KinMixViscosity = ScaleLinker();
	KinMixViscosity:set_import_1(Inverse_RHO)
	

	local EfectiveKinViscosity = nil
	local KinTurbulentViscosity = nil
	local TurbulenceDisc = nil
	if self.turbViscMethod=="no" then
		KinTurbulentViscosity = 0.0
	elseif self.turbViscMethod=="dyn" or self.turbViscMethod=="sma" then
		if self.turbViscMethod == "dyn" then
			KinTurbulentViscosity = FV1DynamicTurbViscData(self.approxSpace,u)
		else
			KinTurbulentViscosity = FV1SmagorinskyTurbViscData(self.approxSpace,u,self.modellconstant)
		end
		KinTurbulentViscosity:set_turbulence_zero_bnd(self.geometry.boundaries)
		KinTurbulentViscosity:set_kinematic_viscosity(0.0)
		
	elseif self.turbViscMethod=="komegaSST" then
		TurbulenceDisc = RANSTurbulenceFV1({"k", "omega"}, Inner_total)
		KinTurbulentViscosity = TurbulenceDisc:turbulent_kinematic_viscosity()
	end



	if self.turbViscMethod=="no" then
		EfectiveKinViscosity = KinMixViscosity
	else
		EfectiveKinViscosity = ScaleAddLinkerNumber()
		EfectiveKinViscosity:add(1.0,KinMixViscosity)
		EfectiveKinViscosity:add(Scale_RHO,KinTurbulentViscosity)

	end




	local MixViscosity = ScaleLinker();
	MixViscosity:set_import_1(Density)
	MixViscosity:set_import_2(EfectiveKinViscosity)



	local DiffusionViscosity = nil
	if self.turbViscMethod=="no" then
		DiffusionViscosity = ConstUserNumber(self.rho_a * self.nu_a)
	else
		DiffusionViscosity = ConstUserNumber(self.rho_a * self.nu_a)
	end


	local TurbulentViscosity = ScaleLinker();
	TurbulentViscosity:set_import_1(AirDensity)
	TurbulentViscosity:set_import_2(KinTurbulentViscosity)

	self.TurbulenceDisc = TurbulenceDisc
	self.TurbulentViscosity = TurbulentViscosity
	self.KinTurbulentViscosity = KinTurbulentViscosity
	self.MixViscosity = MixViscosity
	self.KinMixViscosity = KinMixViscosity
	self.EfectiveKinViscosity = EfectiveKinViscosity
	
	---------------------------------------------------------------------- Normal
	
	local ss_value = math.atan(self.FricMu_2)*180/3.1415926
	
	local Normal = nil
	if self.boolAveNormal then
		Normal = DuneNormal(self.approxSpace,u)
		Normal:set_theta(ss_value)
		Normal:set_gradient_limit(self.grad_limit)
		Normal:set_phase_parameters(InterfaceValues)
	else
		Normal = InterfaceNormalLinker()
	end
	
	self.Normal = Normal

	---------------------------------------------------------------------- Sediment Velocity

	local iter = 0
	local Vs = InterfaceValues:RelVel_ext(0.0,self.rho_a,self.dp,self.rho_s,self.gravity)
	local re = InterfaceValues:RE(self.nu_a*self.rho_a,self.rho_a,self.dp,Vs)
	local Cd = InterfaceValues:CD(re,self.drag_mod)


	--RelVel = ScaleAddLinkerVector()
	--RelVel:add(-Vs, self.Normal)
	local RelVel = RelativeVelocity(self.approxSpace,u)
	RelVel:set_phase_parameters(InterfaceValues)

	
	self.RelVel = RelVel

	---------------------------------------------------------------------- Diffusion (Erotion)

	local Diff_beta = 0.217
	local d1= 250e-06
	local A1= 1.673
	local k0 = 1.0 + A1 * (1-self.dp/d1)
	local Diff_factor = Diff_beta*k0*k0*Vs/(2*self.rho_a*(math.abs(self.gravity)))
	local Diffusion = GranularDiffusionLinker();
	Diffusion:set_mix_viscosity(DiffusionViscosity)
	Diffusion:set_diff_factor(Diff_factor)
	Diffusion:set_phase_parameters(InterfaceValues)
	
	
	--local Dreg = ConstUserMatrix(0.2*Vs)


	self.Diffusion = Diffusion
	--self.Dreg = Dreg
	
	
	---------------------------------------------------------------------- SaltationFlux
	local SaltFlux = SaltationFluxLinker()
	SaltFlux:set_phase_parameters(InterfaceValues)
	SaltFlux:set_saltation_coefficient(5.5)
	SaltFlux:set_threshold_friction_velocity(0.22)
	SaltFlux:set_normal_epsilon(1e-8)
	SaltFlux:set_delta_epsilon(1e-8)
	SaltFlux:set_velocity_epsilon(1e-8)
	SaltFlux:set_direction_mode(2)
	SaltFlux:set_turbulent_kinematic_viscosity(self.KinTurbulentViscosity)
	SaltFlux:set_normal(self.Normal)
	self.SaltFlux = SaltFlux
	
	---------------------------------------------------------------------- AvalancheFlux
	local AvalancheFlux = AvalancheFluxLinker()
	AvalancheFlux:set_phase_parameters(InterfaceValues)
	AvalancheFlux:set_avalanche_coefficient(1.0)
	AvalancheFlux:set_angle_of_repose_degrees(34.0)
	AvalancheFlux:set_normal_epsilon(1.0e-8)
	AvalancheFlux:set_delta_epsilon(1.0e-8)
	AvalancheFlux:set_vertical_epsilon(1.0e-8)
	self.AvalancheFlux = AvalancheFlux
	
	---------------------------------------------------------------------- Vel-Diffusion (Avalanching)
	
	

	local SlipDiff = nil
	local SlipVel = nil
	if self.boolSlipDiff then
		SlipDiff = SlipDiffusion(self.approxSpace,u)
		SlipDiff:set_theta(ss_value)
		SlipDiff:set_diff(self.SlipVelValue)
		SlipDiff:set_slope_limit(self.slope_limit)
		SlipDiff:set_normal(Normal)
		SlipDiff:set_phase_parameters(InterfaceValues)
		
		if self.boolAveDiff then
			SlipDiff:set_diffusion(Diffusion)
		end
	else if self.boolSlipVel then
			SlipVel = SlipVelocity(self.approxSpace,u)
			SlipVel:set_theta(ss_value)
			SlipVel:set_vel(self.SlipVelValue)
			SlipVel:set_slope_limit(self.slope_limit)
			SlipVel:set_normal(Normal)
			SlipVel:set_phase_parameters(InterfaceValues)
		end
	end
	
	self.SlipDiff = SlipDiff
	self.SlipVel = SlipVel
	
	
	print("Clousures Discretization DONE")


end
--------------------------------------------------------------------------------
-- ConnectClosures
--------------------------------------------------------------------------------
myProblem.ConnectClosures = function(self, NavierStokesDisc)
    if not (self.bStokes) and not(self.boolSaltationFlux) then
		self.Density:set_volume_fraction(NavierStokesDisc:volume_fraction())
	end
    self.Diffusion:set_velocity_gradient(NavierStokesDisc:velocity_grad())
    self.KinMixViscosity:set_import_2(NavierStokesDisc:mix_viscosity())
    if not self.boolAveNormal then
		self.Normal:set_volume_grad(NavierStokesDisc:volume_fraction_grad())
	end
    
	self.SaltFlux:set_volume_grad(NavierStokesDisc:volume_fraction_grad())
	self.SaltFlux:set_velocity_gradient(NavierStokesDisc:velocity_grad())
	self.SaltFlux:set_velocity(NavierStokesDisc:velocity())
	
	self.AvalancheFlux:set_volume_grad(NavierStokesDisc:volume_fraction_grad())

end

------------------------------------------------------------
-- UPDATE STATE-DEPENDENT PARAMETERS
------------------------------------------------------------

myProblem.UpdateParameters = function (self)
	if self.turbViscMethod=="dyn" or self.turbViscMethod=="sma" then
		self.KinTurbulentViscosity:update()
	end
	--self.gamma:update()
	self.RelVel:update()
	if self.boolSlipDiff then
		self.SlipDiff:update()
	else if self.boolSlipVel then
			self.SlipVel:update()
		end
	end
	if self.boolAveNormal then
		self.Normal:update()
	end
end
--------------------------------------------------------------------------------
-- Navier-Stokes Discretization
--------------------------------------------------------------------------------
myProblem.Discretization = function (self)
	local Inner_total = self.geometry.innerSubsets
	local NavierStokesDisc = NavierStokesFV1M (self.nv_cmp_tbl, Inner_total)
	NavierStokesDisc:set_exact_jacobian (self.exactJac)
	NavierStokesDisc:set_stokes (self.bStokes)
	NavierStokesDisc:set_laplace (self.bNoLaplace)
	NavierStokesDisc:set_upwind (self.upwind_m)
	NavierStokesDisc:set_upwind_vol(self.upwind_t)
	NavierStokesDisc:set_peclet_blend (self.bPecletBlend)
	NavierStokesDisc:set_stabilization (self.stab, self.diffLength)
	NavierStokesDisc:set_div_correction (self.div_correction)
	NavierStokesDisc:set_transport_ip_velocity(self.boolIPVelocity)
	NavierStokesDisc:set_mass_mean(self.boolDensityMean)
	NavierStokesDisc:set_mass_term(self.boolMassTerm)
	NavierStokesDisc:set_saltation_riemann(true)
	NavierStokesDisc:set_viscosity_averaging("interpolated")
	
	
	NavierStokesDisc:set_kinematic_viscosity (self.EfectiveKinViscosity)
	--NavierStokesDisc:set_average_gamma(self.gamma)
	if self.boolSaltationFlux then
		NavierStokesDisc:set_saltation_flux(self.SaltFlux)
	end
	
	if self.timeMethod == "limex" and self.boolMassTerm and not(self.bStokes) then
		NavierStokesDisc:set_density(self.Density,false)
		NavierStokesDisc:set_limex_correction(true)
	else
		NavierStokesDisc:set_density(self.Density,true)
	end

		
	if self.boolRelativeVel then
		NavierStokesDisc:set_relative_velocity(self.RelVel,self.riemman)
		NavierStokesDisc:set_upwind_rel(self.upwind_r)
	end
	if self.boolSlipDiff then
		NavierStokesDisc:set_diffusion(self.SlipDiff)
	else
		if self.boolAveDiff then
			NavierStokesDisc:set_diffusion(self.Diffusion)
		end
		if self.boolSlipVel then
			NavierStokesDisc:set_slip_velocity(self.SlipVel)
		end
		
	end
	if(self.boolpress_jump) then
		NavierStokesDisc:set_pressure_jump ( self.diffLength)
	end
	
	NavierStokesDisc:set_phase_parameters(self.InterfaceValues)
	
	print("Space Discretization DONE")
	
	self.NavierStokesDisc = NavierStokesDisc
	
	local TurbulenceDisc = nil
	if self.turbViscMethod =="komegaSST" then
		TurbulenceDisc = self:TurbulenceDiscretization()
	end
	
	
	return NavierStokesDisc, TurbulenceDisc
end

--------------------------------------------------------------------------------
-- k-Omega SST Discretization
--------------------------------------------------------------------------------

myProblem.TurbulenceDiscretization = function (self)

	self.TurbulenceDisc:set_velocity(self.NavierStokesDisc:velocity_ip())
	self.TurbulenceDisc:set_upwind("full")
	self.TurbulenceDisc:set_velocity_gradient(self.NavierStokesDisc:velocity_grad())
	self.TurbulenceDisc:set_wall_distance("WallDistance")
	self.TurbulenceDisc:set_kinematic_viscosity(self.KinMixViscosity)

	self.TurbulenceDisc:set_linearize_turbulent_viscosity(true)
	self.TurbulenceDisc:set_linearize_destruction_coupling(true)
	self.TurbulenceDisc:set_linearize_f2(true)
	self.TurbulenceDisc:set_cross_diffusion_linearization(0)
	self.TurbulenceDisc:set_linearize_exported_viscosity(true)
	
	print("Space Discretization (k-OMEGA SST) DONE")
	

	return self.TurbulenceDisc
end
--------------------------------------------------------------------------------
-- BoundaryConditions
--------------------------------------------------------------------------------
myProblem.CreateBoundaryConditions = function(self, NavierStokesDisc, geometry)
    local boundaries = {}
    local inlet = NavierStokesInflowFV1M(NavierStokesDisc)
	local dirichlet = DirichletBoundary()

    if self.problem == "flow" then
        inlet:add("SandInflow", "SandInflow", "Left,Top")
        if self.dim == 3 then
			inlet:add("SandInflow", "SandInflow", "Back,Front")
        end
		table.insert(boundaries, inlet)
		
		local outlet = NavierStokesNoNormalStressOutflowFV1M(NavierStokesDisc)
        outlet:add("Right")
        outlet:set_phase_parameters(self.InterfaceValues)
        table.insert(boundaries, outlet)
        
		local wall = NavierStokesWall(NavierStokesDisc)
		wall:add("Bottom")
		table.insert(boundaries, wall)
		
		dirichlet:add(0.0, "c", "Left,Top")
		if self.dim == 3 then
			dirichlet:add(0.0, "c", "Back,Front")
        end
		
        table.insert(boundaries, dirichlet)
	else
		inlet:add("SandInflow", "SandInflow", geometry.boundaries)
		table.insert(boundaries, inlet)
		
		local neumann = NeumannBoundaryFV1("c")
        neumann:add("SandTopFlux", "Top", "Inner")
        neumann:add(0.0, "Bottom", "Inner")

        if self.dim == 2 then
            dirichlet:add(0.0, "c", "Left")

            if self.simCaseBnd == 1 then
                neumann:add(0.0, "Right", "Inner")
            elseif self.simCaseBnd == 2 then
                dirichlet:add(0.0, "c", "Right")
            else
                error("Unknown avalanche boundary case")
            end
        else
            dirichlet:add(1.0, "c", "Bottom")

            if self.simCaseBnd == 1 then
                dirichlet:add(0.0, "c", "Left1,Left2,Front1,Front2,Right,Back1,Back2,Back3")
            elseif self.simCaseBnd == 2 then
                dirichlet:add(0.0, "c", "Left2,Front1,Front2,Right,Back1,Back2,Back3")
                neumann:add(0.0, "Left1,Back1", "Inner")
            else
                error("Unknown avalanche boundary case")
            end
        end

        table.insert(boundaries, dirichlet)
        table.insert(boundaries, neumann)
	end


    if self.turbViscMethod=="komegaSST" then
		local RANSBnd = DirichletBoundary()
		if self.problem == "flow" then
			RANSBnd:add("InletK","k", "Left,Top")
			RANSBnd:add("InletOmega", "omega", "Left,Top")
			if self.dim == 3 then
				RANSBnd:add("InletK","k", "Back,Front")
				RANSBnd:add("InletOmega", "omega", "Back,Front")
			end
			
			RANSBnd:add(0.0, "k", "Bottom")
			RANSBnd:add(self.omegaWall, "omega", "Bottom")
			
			RANSBnd:add(self.kIn, "k", "Right")
			RANSBnd:add(self.omegaIn, "omega", "Right")
			
			
		else
			RANSBnd:add(0.0, "k", geometry.boundaries)
			RANSBnd:add(0.0, "omega", geometry.boundaries)
		end

		table.insert(boundaries, RANSBnd)
	end

    return boundaries
end

--------------------------------------------------------------------------------
-- Time Discretization
--------------------------------------------------------------------------------
myProblem.TimeDiscretization = function (self,domainDisc)
	
	print("Time Discretization")
	-- create the assembled operator for the solver
	local timeDisc = nil
	-- create time discretization
	if self.timeMethod=="cn" then
		timeDisc = ThetaTimeStep(domainDisc)
		timeDisc:set_theta(0.5) -- Crank-Nicolson method
	end
	if self.timeMethod=="euler" then
		timeDisc = ThetaTimeStep(domainDisc)
		timeDisc:set_theta(1) -- implicit Euler
	end
	if self.timeMethod=="fracstep" then
		timeDisc = ThetaTimeStep(domainDisc,"FracStep")
	end
	if self.timeMethod=="alex" then
		timeDisc = ThetaTimeStep(domainDisc, "Alexander")
	end
	
	print("Time Discretization DONE")
	
	self.timeDisc = timeDisc
	return timeDisc
end

--------------------------------------------------------------------------------
-- SOLVER
--------------------------------------------------------------------------------
myProblem.Smoother = function (self,smoother_name,ordering_name)
	--------------
	-- Smoothers
	--------------
	-- base solver
	local smoother = nil
	if smoother_name == "ilut" then
		local ilut = ILUT(self.eps_ilut)
		ilut:set_damp(self.damping_ilu)
		ilut:set_sort(false)
		smoother = ilut
	elseif smoother_name == "ilu" then
		local ilu = ILU()
		ilu:set_beta(self.value_beta)
		ilu:set_damp(self.damping_ilu)

		if ordering_name == "native" then
			ilu:set_ordering_algorithm(NativeCuthillMcKeeOrdering())
		elseif ordering_name == "boost" then
			ilu:set_ordering_algorithm(BoostCuthillMcKeeOrdering())
		elseif ordering_name == "minimum" then
			ilu:set_ordering_algorithm(BoostMinimumDegreeOrdering())
		elseif ordering_name == "scc" then
			ilu:set_ordering_algorithm(SCCOrdering())
		elseif ordering_name == "none" then
			ilu:set_sort(false)
		else
			error("Unknown ILU ordering: " .. tostring(ordering_name))
		end

		ilu:set_inversion_eps(1.e-12)
		ilu:enable_consistent_interfaces(true)
		ilu:enable_overlap(false)
		smoother = ilu
	elseif smoother_name == "jac" then
		local jac = Jacobi (0.7);
		smoother = jac
	elseif smoother_name == "bgs" then
		local bgs = BlockGaussSeidel ();
		smoother = bgs
	elseif smoother_name == "gs" then
		local gs = GaussSeidel()
		gs:enable_consistent_interfaces(false)
		gs:enable_overlap(true)
		smoother = gs
	elseif smoother_name == "sgs" then
		local sgs = SymmetricGaussSeidel ()
		sgs:enable_consistent_interfaces(true)
		sgs:enable_overlap(true)
		smoother = sgs
	elseif smoother_name == "egs" then
		local egs = ElementGaussSeidel();
		smoother = egs
	elseif smoother_name == "cgs" then
		local cgs = ComponentGaussSeidel(1.0, {"p"})
		cgs:set_alpha(1.0)
		cgs:set_beta(1.0)
		cgs:set_weights(true)
		smoother = cgs
	else
		print ("Smoother not defined"); exit();
	end
	return smoother
end
myProblem.BaseSolver = function (self,baseSolver_name)
	--------------
	-- Base Solver
	--------------
	local Solver
	if baseSolver_name == "SuperLu" then
		--local Solver = LU()
		Solver = AgglomeratingSolver(SuperLU());
	else
		local baseConvCheck=ConvCheck(1000, 1e-8, 1e-03, false)
		Solver = self:LinearSolver("bicgStab",baseConvCheck,"ilu")
	end
	return Solver
end
	

myProblem.Preconditioner = function (self,precond_name,smoother_name, ordering_name,pre_smooth_num, post_smooth_num, baseSolverG, baseLevel)
	print("Creating Precond")
	local preconditioner = nil
	local smoother = self:Smoother(smoother_name,ordering_name)
	if precond_name == "gmg" then
		local baseSolver = self:BaseSolver(baseSolverG)
		local gmg = GeometricMultiGrid(self.approxSpace)
		gmg:set_discretization(self.domainDisc)
		gmg:set_base_level(baseLevel)
		gmg:set_base_solver(baseSolver)
		gmg:set_smoother(smoother)
		gmg:set_cycle_type(1)
		gmg:set_num_presmooth(pre_smooth_num)
		gmg:set_num_postsmooth(post_smooth_num)
		gmg:set_rap( self.rap)
		gmg:set_smooth_on_surface_rim(self.rim)
		
		local mgStats = MGStats()
		mgStats:set_filename_prefix("GMG_debug")
		mgStats:set_exit_on_error(false)
		mgStats:set_write_err_vecs(false)
		mgStats:set_write_err_diffs(false)
		gmg:set_mg_stats(mgStats)

		-- gmg:set_damp(MinimalResiduumDamping())
		gmg:set_damp(self.damping_gmg)
		-- gmg:set_damp(MinimalEnergyDamping())
		--transfer = StdTransfer()
		--transfer:set_restriction_damping(1.0)
		--transfer:set_prolongation_damping(0.6)
		
		--gmg:set_transfer(transfer)
		preconditioner = gmg
	else
		preconditioner = smoother
	
	end
	
	return preconditioner
end

myProblem.LinearSolver = function (self,LinSolver_name,LinearConvCheck,precond_name,smoother_name, ordering_name, pre_smooth_num, post_smooth_num, baseSolverG,baseLevel)
	local LinearSolver = nil
	local preconditioner = self:Preconditioner(precond_name,smoother_name, ordering_name,pre_smooth_num, post_smooth_num,baseSolverG,baseLevel)

	if LinSolver_name == "gmres" then
		-- create Linear Solver
		LinearSolver = GMRES(10)
		LinearSolver:set_preconditioner(preconditioner)
	elseif LinSolver_name == "bicgStab" then
		LinearSolver = BiCGStab()
		LinearSolver:set_preconditioner(preconditioner)
	elseif LinSolver_name == "gmg" then
		LinearSolver = LinearSolver()
		LinearSolver:set_preconditioner(preconditioner)
	elseif LinSolver_name == "ilu" then
		LinearSolver = LinearSolver()
		LinearSolver:set_preconditioner(preconditioner)
	else
		print ("LinearSolver not defined"); exit();
	end
	
	LinearSolver:set_convergence_check(LinearConvCheck)
	
	return LinearSolver
end
myProblem.CreateSolver = function (self, domainDisc)

	self.domainDisc = domainDisc
	local op = nil
	
	--------------
	-- LineSearch
	--------------
	local NewtonLineSearch = nil
	if true then
		NewtonLineSearch = StandardLineSearch()
		NewtonLineSearch:set_maximum_steps(self.lambdamaxSteps)
		NewtonLineSearch:set_lambda_start(self.lambdaStart)
		NewtonLineSearch:set_reduce_factor(0.5)
		NewtonLineSearch:set_accept_best(true)
		NewtonLineSearch:set_check_all(false)
		NewtonLineSearch:set_suff_descent_factor(0.2)
		NewtonLineSearch:set_maximum_defect(2e20)
	else
		NewtonLineSearch = TrustRegionMethod()
		NewtonLineSearch:set_maximum_steps(3)
		NewtonLineSearch:set_lambda_start(self.lambdaStart)
		NewtonLineSearch:set_reduce_factor(0.5)
		NewtonLineSearch:set_accept_best(true)
		NewtonLineSearch:set_check_all(false)
		NewtonLineSearch:set_suff_descent_factor(0.25)
		NewtonLineSearch:set_maximum_defect(2e20)
	end
	----------------------
	-- NoLinear COnvCheck
	----------------------
	
	local NewtonSteadyConvCheck = CompositeConvCheck(self.approxSpace)
	NewtonSteadyConvCheck:set_maximum_steps(self.max_newton_steps_steady_state)
	NewtonSteadyConvCheck:set_group_check(self.vel_cmp_tbl, self.SteadyAbsDefect, self.SteadyRedDefect)
	NewtonSteadyConvCheck:set_component_check("p", self.SteadyAbsDefect, self.SteadyRedDefect)
	NewtonSteadyConvCheck:set_component_check("c", self.SteadyAbsDefect, self.SteadyRedDefect)
	if self.turbViscMethod == "komegaSST" then
		NewtonSteadyConvCheck:set_component_check("k", self.SteadyAbsDefect, self.SteadyRedDefect)
		NewtonSteadyConvCheck:set_component_check("omega", self.SteadyAbsDefect, self.SteadyRedDefect)
	end
	NewtonSteadyConvCheck:disable_rest_check()
	NewtonSteadyConvCheck:set_verbose(true)
	
	local NewtonConvCheck = CompositeConvCheck(self.approxSpace)
	NewtonConvCheck:set_maximum_steps(self.max_newton_steps_transient)
	NewtonConvCheck:set_group_check(self.vel_cmp_tbl, self.AbsDefect, self.RedDefect)
	NewtonConvCheck:set_component_check("p", self.AbsDefect, self.RedDefect)
	NewtonConvCheck:set_component_check("c", self.AbsDefect, self.RedDefect)
	if self.turbViscMethod == "komegaSST" then
		NewtonConvCheck:set_component_check("k", self.AbsDefect, self.RedDefect)
		NewtonConvCheck:set_component_check("omega", self.AbsDefect, self.RedDefect)
	end
	NewtonConvCheck:disable_rest_check()
	NewtonConvCheck:set_verbose(true)
	
	local LinearConvCheckImp=ConvCheck(self.max_linear_steps_Imp, self.LinAbsDefectImp, self.LinRedDefectImp, true)
	local LinearConvCheckLim=ConvCheck(self.max_linear_steps_Lim, self.LinAbsDefectLim, self.LinRedDefectLim, true)
	local LimexConvCheck=ConvCheck(1, 1e-12, 1e-12, true)
	LimexConvCheck:set_supress_unsuccessful(true)
	

	
	-----------------
	-- Linear Solver
	-----------------
	print("Creating linear solver")
	local LinearSolverLim = self:LinearSolver(self.linear_solver_Lim,LinearConvCheckLim,self.precondLim,self.smootherLim, self.orderingLim, self.pre_smoothLim, self.post_smoothLim, self.baseSolverLim,self.gmgBaseLevelLim)
	local LinearSolverImp = self:LinearSolver(self.linear_solver_Imp,LinearConvCheckImp,self.precondImp,self.smootherImp, self.orderingImp, self.pre_smoothImp, self.post_smoothImp, self.baseSolverImp,self.gmgBaseLevelImp)
	

	

	local NewtonSolverSteady = nil
	if self.doSteadyState then
		NewtonSolverSteady = NewtonSolver()
		NewtonSolverSteady:set_linear_solver(LinearSolverImp)
		NewtonSolverSteady:set_convergence_check(NewtonSteadyConvCheck)
		NewtonSolverSteady:set_line_search(NewtonLineSearch)
		
		if self.NewtonSteadyDebug then
			local dbgWriter_steady = GridFunctionDebugWriter(self.approxSpace)
			dbgWriter_steady:set_vtk_output(true)
			dbgWriter_steady:set_conn_viewer_output(false)
			dbgWriter_steady:set_base_dir(self.debug_dir)
			NewtonSolverSteady:set_debug(dbgWriter_steady)
		end
		
	end
		
	local limex = nil
	local NLSolver = NewtonSolver()
	
	
	if self.timeMethod == "limex" then
		NLSolver:set_linear_solver(LinearSolverLim)
		NLSolver:set_convergence_check(LimexConvCheck)
		NLSolver:auto_update(false)
		limex = self:LimexObject( domainDisc, NLSolver)
	else
		NLSolver:set_linear_solver(LinearSolverImp)
		NLSolver:set_convergence_check(NewtonConvCheck)
		NLSolver:set_line_search(NewtonLineSearch)
		NLSolver:set_reassemble_J_freq(0)
		NLSolver:auto_update(true)
		op = AssembledOperator(self.timeDisc)
		op:init()
		NLSolver:init(op)
		if NLSolver:prepare(self.u) == false then
			print ("Newton solver prepare failed.") return op, NLSolver, NewtonSolverSteady, limex, 0
		end
	end
	
	if self.NewtonDebug then
		local dbgWriter = GridFunctionDebugWriter(self.approxSpace)
		dbgWriter:set_vtk_output(true)
		dbgWriter:set_conn_viewer_output(false)
		dbgWriter:set_base_dir(self.debug_dir)
		NLSolver:set_debug(dbgWriter)
					
	end
	TransientNewtonUpdaterSteady = nil
	if (self.NewtonUpdaterSteady and self.turbViscMethod == "komegaSST")then
		TransientNewtonUpdaterSteady = NewtonUpdaterProjection()
		TransientNewtonUpdaterSteady:add_projection(self.dim+2, 1e-38, 1e20) -- k
		TransientNewtonUpdaterSteady:add_projection(self.dim+3, 1e-10, 1e20) -- omega
		if self.algebraBlockSize == 1 and false then
			local projectionMask = GridFunction(self.approxSpace)
			projectionMask:set(0.0)
			Interpolate(1.0, projectionMask, "c")
			TransientNewtonUpdater:set_projection_mask(projectionMask)
		end
		NewtonSolverSteady:setNewtonUpdater(TransientNewtonUpdaterSteady)
	end
	TransientNewtonUpdater = nil
	if (self.NewtonUpdater and self.timeMethod == "euler")then
	
		TransientNewtonUpdater = NewtonUpdaterProjection()
		TransientNewtonUpdater:set_projection_fct(self.dim+1)
		TransientNewtonUpdater:set_max_threshold(1.1)
		TransientNewtonUpdater:set_min_threshold(-1e-5)
		if self.algebraBlockSize == 1 then
			local projectionMask = GridFunction(self.approxSpace)
			projectionMask:set(0.0)
			Interpolate(1.0, projectionMask, "c")
			TransientNewtonUpdater:set_projection_mask(projectionMask)
		end
		NLSolver:setNewtonUpdater(TransientNewtonUpdater)
	end
	
	
	self.NewtonSolverDescSteady = NewtonSolverDescSteady
	self.NewtonSolverDesc = NewtonSolverDesc
	self.boolSolverDesc = true
	
	print("Solver Setting DONE")
	return op, NLSolver, NewtonSolverSteady, limex
end
--------------------------------------------------------------------------------
-- InitializeSolution
--------------------------------------------------------------------------------
myProblem.InitializeSolution = function(self, u, folder)
    local time = 0.0
    local step = 0
    local time_work_total = 0.0
    local interpolate = true

    if self.boolLoadCheckPoint then
        time, step, time_work_total, interpolate = self:LoadCheckPoint(u, folder)
    end
    
    self.loadedCheckPoint = not interpolate

    if interpolate then
        Interpolate(0.0001, u, "u")
        Interpolate(0.0, u, "v")
        if self.dim == 3 then
			Interpolate(0.0, u, "w")
		end
        Interpolate(0.0, u, "p")
        Interpolate("SandInitialC", u, "c")
		if self.turbViscMethod=="komegaSST" then
			Interpolate("SandInitialK", u, "k")
			Interpolate("SandInitialOmega", u, "omega")
		end
    end

    return time, step, time_work_total
end

--------------------------------------------------------------------------------
-- OutputParameters
--------------------------------------------------------------------------------
myProblem.OutputParameters = function (self)

	local out = VTKOutput()
	out:clear_selection()
	out:select_all(false)
	if self.dim == 2 then
		out:select_nodal ("u,v", "velocity")
	else
		out:select_nodal ("u,v,w", "velocity")
	end
	out:select_nodal ("u", "u")
	out:select_nodal ("v", "v")
	if self.dim == 3 then
		out:select_nodal ("w", "w")
	end
	out:select_nodal ("p", "p")
	out:select_nodal ("c", "c")
	if self.turbViscMethod=="komegaSST" then
		out:select_nodal ("k", "k")
		out:select_nodal ("omega", "omega")
	end
	out:select(self.Density, "Rho")
	out:select(self.NavierStokesDisc:mix_density(), "Rho2")
	out:select(self.NavierStokesDisc:einstein_viscosity(), "Mu_eins")
	out:select(self.NavierStokesDisc:mix_viscosity(), "Mu_I")
	out:select(self.MixViscosity, "MixViscosity")
	out:select(self.TurbulentViscosity, "Mu_t")
	
	out:select(self.RelVel, "RelVel")
	out:select(self.NavierStokesDisc:particle_pressure(), "Ps")
	out:select(self.NavierStokesDisc:particle_pressure_grad(), "DPs")
	out:select(self.NavierStokesDisc:volume_fraction_grad(), "DC")
	out:select(self.NavierStokesDisc:velocity_grad(), "Gamma")
	--out:select(self.gamma, "MeanGamma")
	if (self.boolSlipDiff) then
		out:select_element(self.SlipDiff, "SlipDiff")
	else
		if self.boolSlipVel then
			out:select_element(self.SlipVel, "SVel")
		end
	end
	out:select_element(self.Diffusion, "D")
	out:select_element(self.Normal, "n")
	out:select_element(self.SaltFlux,"SaltFlux" )
	out:select_element(self.AvalancheFlux,"AvalancheFlux" )

	print("Output file setting DONE")
	return out
end
--------------------------------------------------------------------------------
-- CheckPoint
--------------------------------------------------------------------------------
myProblem.SaveCheckPoint = function (self,u,folder)

	if (not(self.boolSaveCheckPoint)) then return end
	
	SaveToFile(u, folder .. "/CheckPoint" .. ".vec")
	print("Saving CheckPoint: DONE")

end

myProblem.LoadCheckPoint = function (self,u,folder)

	
	local boolInterpolate = false
	local filename = folder .. "/Integral.txt"
	local ChekPointname = folder .. "/CheckPoint.vec"
	
	local step = 0
	local time_work_total = 0.0
	local Checkfile = io.open(ChekPointname, "r")

	if Checkfile then
		print("Loading CheckPoint")
		Checkfile:close()
		ReadFromFile(u, ChekPointname)
		
		local file = io.open(filename, "r")
		if file then
			for line in file:lines() do
				-- Skip the header
				if not line:match("^Step") and not line:match("^%-") then

					local columns = {}

					for value in line:gmatch("%S+") do
						table.insert(columns, value)
					end

					-- First column = Step
					step = tonumber(columns[1])

					-- Sixth column = TNSteps
					time_work_total = tonumber(columns[6])
				end
			end
			file:close()
		end
		
		
	else
		print("File does not exist: " .. ChekPointname)
		boolInterpolate = true
	end

	local time = self.DT * step
	
	return time, step, time_work_total, boolInterpolate
end

--------------------------------------------------------------------------------
-- Writing Output parameters
--------------------------------------------------------------------------------

myProblem.WriteValues = function (self, folder, step, time, Value_inner1, Value_inner2, WorkTime, TotalWorkTime, Newton_Steps, Newton_Steps_fail,linsolver_calls,linsolver_steps,boolTotal)

	if(step == 0) then
		local file = io.open(folder .. "/Integral.txt", "w+")
		file:write(string.format("Step\tTime\t\tVol-Dom_1\tVol-Dom_2\tWork time\tTotal work time\tTNSteps\tSNSteps\tFNSteps\tLinCalls LinSteps\n"))
		file:write(string.format("%d\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%d\t%d\t%d\t%d\t%d\n", step, time, Value_inner1, Value_inner2, WorkTime, TotalWorkTime, Newton_Steps, Newton_Steps-Newton_Steps_fail, Newton_Steps_fail,linsolver_calls,linsolver_steps))
		file:close()
		print("Integral file created")
	else
		local file = io.open(folder .. "/Integral.txt", "a")
		if(file) then
			if(boolTotal) then
				file:write(string.format("-----------------------------------------------------------------------------------------------------------\n"))
			end
			file:write(string.format("%d\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%d\t%d\t%d\t%d\t%d\n", step, time, Value_inner1, Value_inner2, WorkTime, TotalWorkTime, Newton_Steps, Newton_Steps-Newton_Steps_fail, Newton_Steps_fail,linsolver_calls,linsolver_steps))
			if(boolTotal) then
				file:write(string.format("-----------------------------------------------------------------------------------------------------------\n"))
			end
			file:close()
		end
	end
end

--------------------------------------------------------------------------------
-- Integrals
--------------------------------------------------------------------------------
myProblem.ComputeIntegrals = function(self, u, innerSubsets)
    local field = self.NavierStokesDisc:volume_fraction()
    local first = Integral(field, u, innerSubsets[1], 0.0)
    local second = 0.0
    if innerSubsets[2] then second = Integral(field, u, innerSubsets[2], 0.0) end
    return first, second
end

--------------------------------------------------------------------------------
-- SteadyState Solution
--------------------------------------------------------------------------------
myProblem.ComputeNonLinearSteadyStateSolution = function(self, u, domainDisc, solver)
	local boolSolution = 1
	-- Fix the mass fraction and solve the linear problem for the momentum
	fixer = DirichletBoundary()
	domainDisc:add(fixer)
	fixer:invert_subset_selection()
	fixer:add("c", "")

	solver:init(AssembledOperator(domainDisc))
	
	solver:prepare(u)

	-- apply the solver for the stationary pressure problem
     print("++++++ STEADY STATE CALCULATION BEGIN ++++++")
	local tBefore_s= os.clock()
	if not solver:apply(u) then
		print("===> THE PREPARATION PHASE FAILED! <===")
		boolSolution = 0
	end
	local tAfter_s = os.clock()
    num_newton_steps = solver:num_newton_steps()
    linsolver_calls = solver:total_linsolver_calls()
    linsolver_steps = solver:total_linsolver_steps()
    average_linear_steps = solver:total_average_linear_steps()
    average_non_linear_rates = solver:total_average_non_linear_rates()
    
    print("num_newton_steps = " .. num_newton_steps .. ".")
    print("linsolver_calls = " .. linsolver_calls .. ".")
    print("linsolver_steps = " .. linsolver_steps .. ".")
    print("average_linear_steps = " .. average_linear_steps .. ".")
    print("average_non_linear_rates = " .. average_non_linear_rates .. ".")
    
    solver:clear_average_convergence();
    
	time_work_steady = tAfter_s-tBefore_s
	print("Computation for steady state took " .. time_work_steady .. " seconds.")
	domainDisc:remove (fixer)
	print("++++++++++++++++++++++++ INITIAL CONDITIONS  (STEADY STATE DONE) ++++++++++++++++++++++++")
	
	if(self.max_newton_steps_steady_state == num_newton_steps) then boolSolution = 1 end
	
	return time_work_steady, linsolver_calls, linsolver_steps, boolSolution

end

--------------------------------------------------------------------------------
-- Solution of NonLinear Problem  (Euler temporal discretization)
--------------------------------------------------------------------------------
myProblem.SolveNonlinearProblem = function (self, u, solver, op, solTimeSeries, dt, step, StartTime, EndTime)

	                    
	local red_factor_fail = self.red_factor_fail
	local red_factor_success = self.red_factor_success
	local incr_factor = self.incr_factor
	local optimal_newton_steps = self.optimal_newton_steps
	local DTmin = self.DTmin
	local DTmax = self.DTmax
	local modifyDT = self.modifyDT
	local StepDebug = self.StepDebug
	local AbsDefect = self.AbsDefect
	local RedDefect = self.RedDefect
	local maxConvRate = self.maxConvRate
	local minConvRate = self.minConvRate

	if(self.boolSolverDesc == false) then print("SolverDesc Not initialized ") return 0, 0, 0, 0, 0, 0 end
	local solverDesc = self.NewtonSolverDesc
		
    local Newton_Steps = 0
    local Newton_Steps_fail = 0
    local linsolver_calls_step = 0
    local linsolver_steps_step = 0
    local CompletedStep = false
    local time = StartTime
    local time2 = StartTime
    local dt_in = dt
	if(modifyDT == false) then dt_in = DTmax end
        

        
	while CompletedStep==false  do
        Newton_Steps = Newton_Steps+1
		-- choose time step

		do_dt = math.min(dt_in,math.max((time+DTmax-time2), 0.0))
		print("Size of timestep dt: " .. do_dt)
		-- setup time Disc for old solutions and timestep
		self.timeDisc:prepare_step(solTimeSeries, do_dt)
	
		-- prepare newton solver
		if solver:prepare(u) == false then
			print ("Newton solver failed at step "..step.."."); return 0, 0, 0, 0, 0, 0;
		end
	
		-- apply newton solver
            
        print("++++++ TIMESTEP " .. step-1 + (time2+do_dt-time)/DTmax .. " BEGIN ++++++")
		if solver:apply(u)  == false then
            Newton_Steps_fail = Newton_Steps_fail+1
            dt_in = math.max(dt_in*red_factor_fail,0.99999*DTmin)
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   Reducing Timestep in step  " .. step-1 + (time2+do_dt-time)/DTmax .. ".")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   new DT          =    " .. dt_in .. "     Time = " .. time2 .."")
			print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   old DT    		=    " .. do_dt .."")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   DTmax       =    " .. DTmax .. "     DTmin = " .. DTmin .."")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   Rho         =    " .. solver:total_average_non_linear_rates() .. ".")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   Dt_factor   =    " .. red_factor_fail .. ".")
            
			
			if dt_in < DTmin  or modifyDT== false then
				print ("xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx            Time step below minimum. Aborting. Failed at step  " .. step-1+ (time2+do_dt-time)/DTmax .. ".");
				if StepDebug then
				
					convCheck =
						{
							type		= "standard",
							iterations	= 5,
							absolute	= self.AbsDefect,
							reduction	= self.RedDefect,
							verbose		= true
						}

					solver:set_convergence_check(ConvCheck(5, self.AbsDefect, self.RedDefect, true))
					local dbgWriter = GridFunctionDebugWriter(self.approxSpace)
					dbgWriter:set_vtk_output(true)
					dbgWriter:set_base_dir(self.debug_dir .."Steady")
					solver:set_debug(dbgWriter)
					solver:init(op)
					self.timeDisc:prepare_step(solTimeSeries, do_dt)
					if solver:prepare(u) == false then
						print ("Newton solver failed at DEBUG step "..step..".");
					end
					print("++++++ DEBUG STEP  BEGIN ++++++")
					solver:apply(u)
					print ("Newton solver failed at DEBUG step "..step..".");
				end
				
				return 0, 0, 0, 0, 0, 0;

			else
				VecScaleAssign(u, 1.0, solTimeSeries:latest())
			
			end
		else
            
            if  ((time2 + do_dt + (1e-07) * DTmin -time)/DTmax > 1.0 or modifyDT == false) then
            
                time= self.timeDisc:future_time()                                        -- update new time
                                                                                
                oldestSol = solTimeSeries:oldest()                                  -- get oldest solution
                                                                                    
				VecScaleAssign(oldestSol, 1.0, u)                                             -- copy values into oldest solution (we reuse the memory here)
                                                                                    
                solTimeSeries:push_discard_oldest(oldestSol, time)                  -- push oldest solutions with new values to front, oldest sol pointer is poped from end
                
                CompletedStep = true
                
                
            else
                
                time2 = self.timeDisc:future_time()                                      -- update new time
                                                                                    
                oldestSol = solTimeSeries:oldest()                                  -- get oldest solution
                                                                                    
				VecScaleAssign(oldestSol, 1.0, u)                                             -- copy values into oldest solution (we reuse the memory here)
                                                                                    
                solTimeSeries:push_discard_oldest(oldestSol, time2)                 -- push oldest solutions with new values to front, oldest sol pointer is poped from end
                
                CompletedStep = false
                
            
            end
            
            average_non_linear_rates = solver:total_average_non_linear_rates()
            num_newton_steps = solver:num_newton_steps()
            
            if CompletedStep then
                frac_step = 0
                Local_Time = time
                Dt_factor = dt_in/do_dt
            else
                frac_step = -1 + (time2-time)/DTmax
                Local_Time = time2
            end
            
            
            --[[if modifyDT then
                if(average_non_linear_rates>maxConvRate and false) then
                    dt_in = math.max(do_dt*red_factor_success,1.00001*DTmin)
                    print ("-------------------------------------------------------------------Time step decrease at Step " .. step-1 + (time2-time)/DTmax .. ", dt =  " .. dt_in .. ". ")
                else if (CompletedStep== false or Dt_factor < 0.98 ) then
						if(average_non_linear_rates<minConvRate or num_newton_steps<=optimal_newton_steps) then
							dt_in=math.min(incr_factor*dt_in,DTmax)
							print ("-------------------------------------------------------------------Time step increased at Step " .. step-1 + (time2-time)/DTmax ..", dt =  " .. dt_in .. ". ")
						end
                    end
                end
            end]]
            
			if modifyDT then
                if(num_newton_steps>30) then
                    dt_in = math.max(do_dt*red_factor_success,1.00001*DTmin)
                    print ("-------------------------------------------------------------------Time step decrease at Step " .. step-1 + (time2-time)/DTmax .. ", dt =  " .. dt_in .. ". ")
                else if (CompletedStep== false or Dt_factor < 0.98 ) then
						if(average_non_linear_rates<minConvRate or num_newton_steps<=optimal_newton_steps) then
							dt_in=math.min(incr_factor*dt_in,DTmax)
							print ("-------------------------------------------------------------------Time step increased at Step " .. step-1 + (time2-time)/DTmax ..", dt =  " .. dt_in .. ". ")
						end
                    end
                end
            end
            
            
            
            CFL=cflNumber(u,do_dt)                                              -- compute CFL number
            print("DT=" .. do_dt .. "");
            print("Time=" .. Local_Time .. "");
            
            print("++++++ TIMESTEP " .. step + frac_step.. "  END ++++++")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   Successful semi Timestep in step  " .. step + frac_step  .. ".")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   new DT    =    " .. dt_in .. "     Time = " .. Local_Time .."")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   old DT    =    " .. do_dt .."")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   DTmax =    " .. DTmax .. "     DTmin = " .. DTmin .."")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   Dt_factor   =    " .. dt_in/do_dt .. ".")
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   Rho   =    " .. average_non_linear_rates .. ".")
                        

		end
  
        
        linsolver_calls_step = linsolver_calls_step + solver:total_linsolver_calls()
        linsolver_steps_step = linsolver_steps_step + solver:total_linsolver_steps()
        solver:clear_average_convergence();
			
	end

  return Newton_Steps, Newton_Steps_fail, linsolver_calls_step, linsolver_steps_step, 1
end

--------------------------------------------------------------------------------
-- LIMEX
--------------------------------------------------------------------------------
myProblem.LimexObject = function ( self, domainDisc, limexNLSolver)


	-- local refObserver = PlotRefOutputObserver("DirichletValue", vtk) -- now obsolete
	local luaObserver = LuaCallbackObserver()

	-- work-around (waiting for implementation of SmartPtr forward to lua...)
	function luaPostProcess(step, time, currdt)
	  print("LUAPostProcess: "..step..","..time..","..currdt)
	  postProcess(luaObserver:get_current_solution(), step, time, currdt)
	  return 0;
	end
	luaObserver:set_callback("luaPostProcess")

	local dtmax = self.DTmax
	local dtmin = self.DTmin
	local dtlimex = self.DTLimex

	--  Euclidean (algebraic) norm
	--local estimator = Norm2Estimator()



	--print (estimator)
	local limexEstimator = CompositeGridFunctionEstimator()
	
	if self.VelErrorNorm == "L2" then
		limexEstimator:add(L2ComponentSpace("u", 2, self.alphaVel))
		limexEstimator:add(L2ComponentSpace("v", 2, self.alphaVel))
	elseif self.VelErrorNorm == "H1" then
		limexEstimator:add(H1SemiComponentSpace("u", 2 , ConstUserMatrix(self.alphaVel)))
		limexEstimator:add(H1SemiComponentSpace("v", 2 , ConstUserMatrix(self.alphaVel)))
	else
		print ("LimexErrorEstimator for velocity not defined"); exit();
	end
	
	if self.PressErrorNorm == "L2" then
		limexEstimator:add(L2ComponentSpace("p", 2, self.alphaPress))
	elseif self.PressErrorNorm == "H1" then
		limexEstimator:add(H1SemiComponentSpace("p", 2, ConstUserMatrix(self.alphaPress) ))
	else
		print ("LimexErrorEstimator for Pressure not defined"); exit();
	end
	
	if self.VolErrorNorm == "L2" then
		limexEstimator:add(L2ComponentSpace("c", 2, self.alphaVol))
	elseif self.VolErrorNorm == "H1" then
		limexEstimator:add(H1SemiComponentSpace("c", 2, self.alphaVol))
	else
		print ("LimexErrorEstimator for VolumeFraction not defined"); exit();
	end
	if self.turbViscMethod=="komegaSST" then
		limexEstimator:add(L2ComponentSpace("k", 2))
		limexEstimator:add(L2ComponentSpace("omega", 2))
	end
	

	
	limexEstimator:use_strict_relative_norms(1)
	print(limexEstimator:config_string())

	-- descriptor for integrator
	local limexDesc = {

	  nstages = self.nstages,
	  steps = {1,2,3,4,5,6},
	  nthreads = 1,
	  domainDisc=domainDisc,
	  nonlinSolver = limexNLSolver,
	  -- makeConsistent = true,
	  
	  tol = self.tol,
	  dt = dtlimex,
	  dtmax = dtmax,
	  dtmin = dtmin,
	  rhoSafetyOPT = 0.25,
	  dtred = self.red_factor_fail,
	  dtIncr = self.incr_factor,
	  matrixCache = false,
	  conservative = true
	  
	}


	-- setup for time integrator
	local limex = util.limex.CreateIntegrator(limexDesc)

	limex:set_time_step(limexDesc.dt)
	limex:set_dt_min(limexDesc.dtmin)
	limex:set_dt_max(limexDesc.dtmax)
	limex:set_reduction_factor(limexDesc.dtred)
	limex:set_increase_factor(limexDesc.dtIncr)
	limex:add_error_estimator(limexEstimator)
	limex:set_stepsize_greedy_order_factor(1.0)
	limex:select_cost_strategy(LimexNonlinearCost())


	if (false) then
		
		--limex:attach_observer(vtkObserver)
		limex:attach_observer(luaObserver)
	end


	--limex:attach_observer(refObserver)

	print ("dtLimex   = "..dtlimex)
	print ("tolLimex  = "..self.tol)
	return limex
end
--------------------------------------------------------------------------------
-- SolutionNonLinearProblem LIMEX
--------------------------------------------------------------------------------
myProblem.SolveNonlinearProblemLimex = function (self, u, limex, NLSolver, time_step, StartTime, EndTime, NewtonLimexSteps)
	
	if(time_step == 1) then
		limex:set_dt_min(1e-010)
	end
	limex:set_start_step(1)
	local fail = limex:apply(u, EndTime, u, StartTime)
	
	print("LIMEX total   = ", limex:get_total_steps())
	print("LIMEX success = ", limex:get_success_steps())
	print("LIMEX fail    = ", limex:get_failed_steps())
	print("LIMEX step    = ", limex:get_step())
	
	n_step = limex:get_step()-1
	
    local Newton_Steps = limex:get_total_steps()
	local Newton_Steps_fail = limex:get_failed_steps()
	linsolver_calls_step = NLSolver:total_linsolver_calls()
	linsolver_steps_step =  NLSolver:total_linsolver_steps()
	limex:set_time_step(self.DTmax/(math.max(n_step-1,1)))
	
	if(time_step == 1) then
		limex:set_dt_min(self.DTmin)
	end
	
	NLSolver:clear_average_convergence();
	if (fail) then
		return Newton_Steps, Newton_Steps_fail, linsolver_calls_step, linsolver_steps_step, 1
	else
		return Newton_Steps, Newton_Steps_fail, linsolver_calls_step, linsolver_steps_step, 0
	end
end

myProblem.CheckJacobian = function(self, domainDisc, u, perturbFct, solveLinearSystem)

	------------------------------------------------------------
	-- Settings
	------------------------------------------------------------

	if solveLinearSystem == nil then
		solveLinearSystem = true
	end

	local linSolver = self.LinearSolver

	local innerSubsets = table.concat(self.geometry.innerSubsets, ",")
	local boundarySubsets = self.geometry.boundaries

	local epsList = {
		1.0e-5,
		1.0e-6,
		1.0e-7,
		1.0e-8,
		1.0e-9,
		1.0e-10
	}


	------------------------------------------------------------
	-- UPDATE STATE-DEPENDENT PARAMETERS
	--
	-- Add here all quantities that must be updated whenever
	-- the current solution changes.
	------------------------------------------------------------

	local function UpdateParameters()
		if params.turbViscMethod=="dyn" or params.turbViscMethod=="sma" then
			self.KinTurbulentViscosity:update()
		end
	end


	------------------------------------------------------------
	-- Jacobian perturbation function
	--
	-- Must be GLOBAL because UG4 Interpolate searches for the
	-- callback by its string name in the global Lua namespace.
	------------------------------------------------------------

	function JacobianPerturbation(x, y, t)
		return x
	end


	------------------------------------------------------------
	-- Build perturbation direction v
	------------------------------------------------------------

	local v = u:clone()
	v:set(0.0)

	Interpolate("JacobianPerturbation", v, perturbFct, innerSubsets)


	------------------------------------------------------------
	-- Working vectors
	------------------------------------------------------------

	local rNewton = u:clone()
	rNewton:set(0.0)

	local r0 = u:clone()
	r0:set(0.0)

	local r1 = u:clone()
	r1:set(0.0)

	local Jv = u:clone()
	Jv:set(0.0)

	local fd = u:clone()
	fd:set(0.0)

	local err = u:clone()
	err:set(0.0)


	------------------------------------------------------------
	-- Adjust current solution
	------------------------------------------------------------

	domainDisc:adjust_solution(u)


	------------------------------------------------------------
	-- Newton defect
	--
	-- Newton computes the defect BEFORE the step update.
	------------------------------------------------------------

	domainDisc:assemble_defect(rNewton, u)


	------------------------------------------------------------
	-- Update state-dependent parameters
	--
	-- This corresponds to the Newton step update before the
	-- Jacobian is assembled.
	------------------------------------------------------------

	UpdateParameters()


	------------------------------------------------------------
	-- Assemble Jacobian
	------------------------------------------------------------

	local J = AssembledLinearOperator(domainDisc)
	domainDisc:assemble_jacobian(J, u)


	------------------------------------------------------------
	-- Assemble the reference defect for the FD test
	--
	-- IMPORTANT:
	-- This is assembled AFTER UpdateParameters(), so that
	--
	--     R(u), J(u), R(u + eps*v)
	--
	-- all use exactly the same frozen parameters.
	------------------------------------------------------------

	r0:set(0.0)
	domainDisc:assemble_defect(r0, u)


	------------------------------------------------------------
	-- Component spaces
	------------------------------------------------------------

	local spaceUV = GridFunctionComponentSpace("u,v")
	local spaceP = GridFunctionComponentSpace("p")
	local spaceC = GridFunctionComponentSpace("c")
	local spaceK = GridFunctionComponentSpace("k")
	local spaceOmega = GridFunctionComponentSpace("omega")

	local spaceKInner = GridFunctionComponentSpace("k", innerSubsets)
	local spaceOmegaInner = GridFunctionComponentSpace("omega", innerSubsets)

	local spacePertInner = GridFunctionComponentSpace(perturbFct, innerSubsets)


	------------------------------------------------------------
	-- Optional Newton / linear-system analysis
	------------------------------------------------------------

	if solveLinearSystem then

		--------------------------------------------------------
		-- Actual Newton correction
		--
		-- J * corr = R
		-- u_new = u - lambda * corr
		--------------------------------------------------------

		local corr = u:clone()
		corr:set(0.0)

		linSolver:init(J, corr)

		if not linSolver:apply(corr, rNewton) then

			print("")
			print("================================================")
			print(" LINEAR SOLVE FOR NEWTON CORRECTION FAILED")
			print("================================================")
			print("")

		else

			print("")
			print("================================================")
			print(" ACTUAL NEWTON CORRECTION")
			print("================================================")

			print(string.format("|corr|         = %.6e", VecNorm(corr)))
			print(string.format("|corr u,v|     = %.6e", spaceUV:norm(corr)))
			print(string.format("|corr p|       = %.6e", spaceP:norm(corr)))
			print(string.format("|corr c|       = %.6e", spaceC:norm(corr)))
			print(string.format("|corr k|       = %.6e", spaceK:norm(corr)))
			print(string.format("|corr omega|   = %.6e", spaceOmega:norm(corr)))
			print(string.format("|corr k inner| = %.6e", spaceKInner:norm(corr)))
			print(string.format("|corr w inner| = %.6e", spaceOmegaInner:norm(corr)))

			print("================================================")
			print("")


			----------------------------------------------------
			-- Norm calls may change parallel storage.
			----------------------------------------------------

			corr:enforce_consistent_type()
			u:enforce_consistent_type()


			----------------------------------------------------
			-- Save original solution
			----------------------------------------------------

			local uSave = u:clone()

			corr:enforce_consistent_type()
			uSave:enforce_consistent_type()


			----------------------------------------------------
			-- Newton line-search diagnostic
			----------------------------------------------------

			local lambdaList = {
				1.0,
				0.5,
				0.25,
				0.125,
				0.0625,
				0.03125,
				0.015625,
				0.0078125,
				0.00390625,
				0.001953125,
				0.0009765625,
				0.00048828125,
				0.000244140625,
				0.0001220703125
			}

			print("NEWTON TRIAL DEFECTS")

			for _, lambda in ipairs(lambdaList) do

				corr:enforce_consistent_type()
				uSave:enforce_consistent_type()
				u:enforce_consistent_type()

				VecScaleAdd2(u, 1.0, uSave, -lambda, corr)

				domainDisc:adjust_solution(u)

				------------------------------------------------
				-- DO NOT call UpdateParameters() here.
				--
				-- The Newton line search uses the parameters
				-- frozen at the step update.
				------------------------------------------------

				local rTrial = u:clone()
				rTrial:set(0.0)

				domainDisc:assemble_defect(rTrial, u)

				print(string.format(
					"lambda = %.6f   |R(u-lambda*corr)| = %.6e",
					lambda,
					VecNorm(rTrial)
				))

				u:enforce_consistent_type()
				uSave:enforce_consistent_type()

				VecAssign(u, uSave)
			end


			----------------------------------------------------
			-- Restore original solution and parameter state
			----------------------------------------------------

			u:enforce_consistent_type()
			uSave:enforce_consistent_type()

			VecAssign(u, uSave)

			domainDisc:adjust_solution(u)

			UpdateParameters()

			print("")
		end

	else

		print("")
		print("================================================")
		print(" LINEAR SYSTEM / NEWTON ANALYSIS DISABLED")
		print("================================================")
		print("")

	end


	------------------------------------------------------------
	-- Re-establish state for Jacobian finite-difference test
	------------------------------------------------------------

	domainDisc:adjust_solution(u)
	UpdateParameters()

	------------------------------------------------------------
	-- Assemble reference defect at the test state
	------------------------------------------------------------

	r0:set(0.0)
	domainDisc:assemble_defect(r0, u)

	------------------------------------------------------------
	-- Assemble a fresh Jacobian at exactly the same test state
	------------------------------------------------------------

	local JTest = AssembledLinearOperator(domainDisc)
	domainDisc:assemble_jacobian(JTest, u)

	------------------------------------------------------------
	-- Compute J*v
	------------------------------------------------------------

	Jv:set(0.0)

	v:enforce_consistent_type()
	Jv:enforce_consistent_type()

	JTest:apply(Jv, v)


	------------------------------------------------------------
	-- Print perturbation information
	------------------------------------------------------------

	print("")
	print("================================================")
	print(" " .. string.upper(perturbFct) .. "-DIRECTION JACOBIAN FINITE-DIFFERENCE TEST")
	print("================================================")
	print("")

	print(string.upper(perturbFct) .. " PERTURBATION v")
	print(string.format("inner : %.6e", spacePertInner:norm(v)))
	print("")


	------------------------------------------------------------
	-- Restore consistent storage after norm calls
	------------------------------------------------------------

	v:enforce_consistent_type()
	u:enforce_consistent_type()
	r0:enforce_consistent_type()
	Jv:enforce_consistent_type()
	
	local rMinus = u:clone()
	rMinus:set(0.0)


	------------------------------------------------------------
	-- Finite-difference Jacobian test
	------------------------------------------------------------

	for _, eps in ipairs(epsList) do

		local uPert = u:clone()

		u:enforce_consistent_type()
		v:enforce_consistent_type()
		uPert:enforce_consistent_type()

		VecScaleAdd2(uPert, 1.0, u, eps, v)

		domainDisc:adjust_solution(uPert)


		--------------------------------------------------------
		-- IMPORTANT:
		--
		-- Do NOT call UpdateParameters() here.
		--
		-- The state-dependent parameters remain frozen at the
		-- state at which J was assembled.
		--------------------------------------------------------

		r1:set(0.0)
		domainDisc:assemble_defect(r1, uPert)
		
		


		--------------------------------------------------------
		local uMinus = u:clone()

		u:enforce_consistent_type()
		v:enforce_consistent_type()
		uMinus:enforce_consistent_type()

		VecScaleAdd2(uMinus, 1.0, u, -eps, v)
		domainDisc:adjust_solution(uMinus)

		rMinus:set(0.0)
		domainDisc:assemble_defect(rMinus, uMinus)
		--------------------------------------------------------


		--------------------------------------------------------
		-- FD = [R(u + eps*v) - R(u)] / eps
		--------------------------------------------------------

		r1:enforce_consistent_type()
		r0:enforce_consistent_type()
		fd:enforce_consistent_type()

		r1:enforce_consistent_type()
		rMinus:enforce_consistent_type()
		fd:enforce_consistent_type()

		VecScaleAdd2(fd, 0.5 / eps, r1, -0.5 / eps, rMinus)


		--------------------------------------------------------
		-- err = FD - J*v
		--------------------------------------------------------

		fd:enforce_consistent_type()
		Jv:enforce_consistent_type()
		err:enforce_consistent_type()

		VecScaleAdd2(err, 1.0, fd, -1.0, Jv)


		--------------------------------------------------------
		-- Global norms
		--------------------------------------------------------

		local normFD = VecNorm(fd)
		local normJv = VecNorm(Jv)
		local normErr = VecNorm(err)

		local denom = math.max(normFD, normJv, 1.0e-30)
		local relErr = normErr / denom
		
		local normFDC = spaceC:norm(fd)
		local normJvC = spaceC:norm(Jv)
		local normErrC = spaceC:norm(err)
		local relErrC = normErrC / math.max(normFDC, normJvC, 1.0e-30)

		print(string.format(
			"C BLOCK: eps = %.1e   |FD| = %.6e   |Jv| = %.6e   |err| = %.6e   rel = %.6e",
			eps, normFDC, normJvC, normErrC, relErrC
		))

		local jacOut = VTKOutput()
		jacOut:clear_selection()
		jacOut:select_all(false)
		jacOut:select_nodal("c", "c")
		
		r1:enforce_consistent_type()
		rMinus:enforce_consistent_type()
		fd:enforce_consistent_type()

		VecScaleAdd2(fd, 0.5 / eps, r1, -0.5 / eps, rMinus)

		fd:enforce_consistent_type()
		Jv:enforce_consistent_type()
		err:enforce_consistent_type()

		VecScaleAdd2(err, 1.0, fd, -1.0, Jv)
		
		if eps == epsList[#epsList] then
			fd:enforce_consistent_type()
			Jv:enforce_consistent_type()
			err:enforce_consistent_type()
			v:enforce_consistent_type()

			jacOut:print_subsets(folder_vtk .. "/Jacobian_FD", fd, self.geometry.allSubsets, 0, 0.0, true)
			jacOut:print_subsets(folder_vtk .. "/Jacobian_Jv", Jv, self.geometry.allSubsets, 0, 0.0, true)
			jacOut:print_subsets(folder_vtk .. "/Jacobian_Error", err, self.geometry.allSubsets, 0, 0.0, true)
			jacOut:print_subsets(folder_vtk .. "/Jacobian_Perturbation", v, self.geometry.allSubsets, 0, 0.0, true)

			print("Jacobian c-block diagnostic files written.")
		end

		--------------------------------------------------------
		-- Detailed block analysis at eps = 1e-6
		--------------------------------------------------------

		if eps == epsList[#epsList] then

			print("")
			print("COMPONENT-WISE JACOBIAN CHECK AT eps = " .. epsList[#epsList])
			print("Input block: " .. perturbFct)
			print("")

			print(string.format(
				"u,v   : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceUV:norm(fd),
				spaceUV:norm(Jv),
				spaceUV:norm(err)
			))

			print(string.format(
				"p     : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceP:norm(fd),
				spaceP:norm(Jv),
				spaceP:norm(err)
			))

			print(string.format(
				"c     : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceC:norm(fd),
				spaceC:norm(Jv),
				spaceC:norm(err)
			))

			print(string.format(
				"k     : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceK:norm(fd),
				spaceK:norm(Jv),
				spaceK:norm(err)
			))

			print(string.format(
				"omega : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceOmega:norm(fd),
				spaceOmega:norm(Jv),
				spaceOmega:norm(err)
			))

			print("")

			print(string.format(
				"k inner     : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceKInner:norm(fd),
				spaceKInner:norm(Jv),
				spaceKInner:norm(err)
			))

			print(string.format(
				"omega inner : |FD| = %.6e   |Jv| = %.6e   |err| = %.6e",
				spaceOmegaInner:norm(fd),
				spaceOmegaInner:norm(Jv),
				spaceOmegaInner:norm(err)
			))

			print("")
		end


		--------------------------------------------------------
		-- Global FD comparison
		--------------------------------------------------------

		print(string.format(
			"eps = %.1e   |FD| = %.6e   |Jv| = %.6e   |FD-Jv| = %.6e   rel = %.6e",
			eps,
			normFD,
			normJv,
			normErr,
			relErr
		))


		--------------------------------------------------------
		-- Restore compatible storage for next epsilon
		--------------------------------------------------------

		u:enforce_consistent_type()
		v:enforce_consistent_type()
		r0:enforce_consistent_type()
		r1:enforce_consistent_type()
		Jv:enforce_consistent_type()
		fd:enforce_consistent_type()
		err:enforce_consistent_type()
	end


	print("================================================")
	print("")

end


return myProblem