From 28279b79e3cc7c278d62532dae2cdea0131b3c2b Mon Sep 17 00:00:00 2001 From: DanGonEsp Date: Mon, 26 May 2025 22:22:47 +0300 Subject: [PATCH] Example for NavierStokes-Transport coupling added --- grids/channel_40x15.ugx | 40 +++ navierstokes-transport_coupling.lua | 461 ++++++++++++++++++++++++++++ 2 files changed, 501 insertions(+) create mode 100644 grids/channel_40x15.ugx create mode 100644 navierstokes-transport_coupling.lua diff --git a/grids/channel_40x15.ugx b/grids/channel_40x15.ugx new file mode 100644 index 0000000..cbc434d --- /dev/null +++ b/grids/channel_40x15.ugx @@ -0,0 +1,40 @@ + + + 0 15.1999999999999993 0 0 0.200000000000000178 0 40 0.200000000000000178 0 40 15.1999999999999993 0 20 15.1999999999999993 0 20 0.200000000000000178 0 0 7.69999999999999929 0 40 7.69999999999999929 0 30 15.1999999999999993 0 10 15.1999999999999993 0 10 0.200000000000000178 0 30 0.200000000000000178 0 20 7.69999999999999929 0 10 7.69999999999999929 0 30 7.69999999999999929 0 0 11.4499999999999993 0 0 3.94999999999999973 0 40 3.94999999999999973 0 40 11.4499999999999993 0 35 15.1999999999999993 0 25 15.1999999999999993 0 15 15.1999999999999993 0 5 15.1999999999999993 0 5 0.200000000000000178 0 15 0.200000000000000178 0 25 0.200000000000000178 0 35 0.200000000000000178 0 20 3.94999999999999973 0 20 11.4499999999999993 0 15 3.94999999999999973 0 5 11.4499999999999993 0 35 3.94999999999999973 0 25 11.4499999999999993 0 5 7.69999999999999929 0 5 3.94999999999999973 0 10 3.94999999999999973 0 15 7.69999999999999929 0 15 11.4499999999999993 0 10 11.4499999999999993 0 25 7.69999999999999929 0 25 3.94999999999999973 0 30 3.94999999999999973 0 35 7.69999999999999929 0 35 11.4499999999999993 0 30 11.4499999999999993 0 + 0 15 15 6 6 16 16 1 2 17 17 7 7 18 18 3 3 19 19 8 8 20 20 4 4 21 21 9 9 22 22 0 1 23 23 10 10 24 24 5 5 25 25 11 11 26 26 2 5 27 27 12 12 28 28 4 5 29 29 13 13 30 30 0 2 31 31 14 14 32 32 4 6 33 33 13 10 34 34 6 13 35 35 10 12 36 36 13 9 37 37 12 13 38 38 9 12 39 39 14 11 40 40 12 14 41 41 11 7 42 42 14 8 43 43 7 14 44 44 8 15 30 33 15 30 33 23 16 34 23 16 34 29 24 35 29 24 35 34 33 35 34 33 35 27 29 36 27 29 36 21 28 37 21 28 37 30 22 38 30 22 38 37 36 38 37 36 38 28 32 39 28 32 39 25 27 40 25 27 40 31 26 41 31 26 41 40 39 41 40 39 41 17 31 42 17 31 42 19 18 43 19 18 43 32 20 44 32 20 44 43 42 44 43 42 44 + 0 15 30 6 33 15 13 30 33 15 33 30 1 23 16 10 34 23 6 16 34 23 34 16 5 29 24 13 35 29 10 24 35 29 35 24 6 34 33 10 35 34 13 33 35 34 35 33 5 27 29 12 36 27 13 29 36 27 36 29 4 21 28 9 37 21 12 28 37 21 37 28 0 30 22 13 38 30 9 22 38 30 38 22 12 37 36 9 38 37 13 36 38 37 38 36 4 28 32 12 39 28 14 32 39 28 39 32 5 25 27 11 40 25 12 27 40 25 40 27 2 31 26 14 41 31 11 26 41 31 41 26 12 40 39 11 41 40 14 39 41 40 41 39 2 17 31 7 42 17 14 31 42 17 42 31 3 19 18 8 43 19 7 18 43 19 43 18 4 32 20 14 44 32 8 20 44 32 44 20 7 43 42 8 44 43 14 42 44 43 44 42 + + + 12 13 14 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 + 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 + 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 + + + 8 9 19 20 21 22 4 + 8 9 10 11 12 13 14 15 + + + 2 3 7 17 18 + 4 5 6 7 + + + 5 10 11 23 24 25 26 + 16 17 18 19 20 21 22 23 + + + 0 1 6 15 16 + 0 1 2 3 + + + + + + + + 7 1 + + + 0 0 + + + diff --git a/navierstokes-transport_coupling.lua b/navierstokes-transport_coupling.lua new file mode 100644 index 0000000..b16b310 --- /dev/null +++ b/navierstokes-transport_coupling.lua @@ -0,0 +1,461 @@ +------------------------------------------------------------------------------------------ +-- Navier-Stokes equation, 2d +-- Discretization: Vertex-centered, stabilized +------------------------------------------------------------------------------------------ + +-- Load utility scripts (e.g. from from ugcore/scripts) +ug_load_script ("ug_util.lua") +ug_load_script ("util/load_balancing_util.lua") + +-- Required Plugins +PluginRequired("NavierStokes") +PluginRequired("ConvectionDiffusion") + + +------------------------------------------------------------------------------------------ +-- Get command line parameters +------------------------------------------------------------------------------------------ + +-- Physical parameters + +nu_a = util.GetParamNumber("-visc_a", 1.48e-05, "kinematic viscosity") +rho_a = util.GetParamNumber("-rho_a", 1.2, "Air Density") +inflow = util.GetParamNumber("-inflow", 1.0, "max. inflow velocity") +bStokes = util.HasParamOption("-Stokes", "If defined, only Stokes Eq. computed") +c_init = util.GetParamNumber("-initial concentration", 1.0, "max volume fraction") + +-- Numerical parameters of the discretization +numRefs = util.GetParamNumber("-numRefs", 4, "number of grid refinements") +numPreRefs = util.GetParamNumber("-numPreRefs", 2, "number of prerefinements (parallel)") +bNoLaplace = util.HasParamOption("-noLaplace", "If defined, only laplace term used") +bExactJac = util.HasParamOption("-exactJac", "If defined, exact jacobian used") +bPecletBlend= util.HasParamOption("-PecletBlend", "If defined, Peclet Blend used") +upwind = util.GetParam("-upwind", "full", "Upwind type (no, full, weighted, lps, pos, reg)") +bPac = util.HasParamOption("-pac", "If defined, pac upwind used") +stab = util.GetParam("-stab", "fields", "Stabilization type (fields or flow)") +diffLength = util.GetParam("-difflength", "cor", "Diffusion length type (raw, fivepoint or cor)") + +dt = util.GetParamNumber("-dt",10, "TimeStep") +timeMethod = util.GetParam("-timeMethod","euler") +numTimeSteps = util.GetParamNumber("-numTimeSteps", 10 ) +outputFactor = util.GetParam("-output", 1, "output every ... steps") + + +-- Parameters of the solver +ilu_beta = util.GetParamNumber("-iluBeta", -0.2, "choose a negative value depending on the convection rate") + + +------------------------------------------------------------------------------------------ +-- Geometry data and export file name +------------------------------------------------------------------------------------------ + +dim = 2 +gridName = "grids/channel_40x15.ugx" +file_name ="Transport" + +-- Subsets used in the problem +allSubsets = "Inner,Left, Right,Top, Bottom," + + +------------------------------------------------------------------------------------------ +-- Print the parameters +------------------------------------------------------------------------------------------ + +print (" Geometry: " .. gridName .. ", dim = " .. dim) +print (" Physical parameter:") +print (" visc = " .. nu_a) +print (" inflow = " .. inflow) +print (" Stokes = " .. tostring (bStokes)) +print (" Numerical parameter of the discretization:") +print (" numRefs = " .. numRefs) +print (" numPreRefs = " .. numPreRefs) +print (" noLaplace = " .. tostring (bNoLaplace)) +print (" exactJac = " .. tostring (bExactJac)) +print (" PecletBlend = " .. tostring (bPecletBlend)) +print (" upwind = " .. upwind) +print (" pac = " .. tostring (bPac)) +print (" stab = " .. stab) +print (" difflength = " .. diffLength) +print (" Numerical parameter of the solvers:") +print (" iluBeta = " .. ilu_beta) + + + +------------------------------------------------------------------------------------------ +-- Initialize UG4, load, refine and distribute the grid +------------------------------------------------------------------------------------------ + +InitUG (dim, AlgebraType("CPU", dim +2)) + +if dim == 3 then + fct_cmp_tbl = {"u", "v", "w", "p"} + vel_cmp_tbl = {"u", "v", "w"} +else + fct_cmp_tbl = {"u", "v", "p"} + vel_cmp_tbl = {"u", "v"} +end + +-- Create the domain, load the grid and refine it +dom = util.CreateDomain (gridName, numPreRefs) +balancer.RefineAndRebalanceDomain (dom, numRefs - numPreRefs) + +print ("Domain-info:") +print (dom:domain_info():to_string()) + +-- Create the vertex-centered approximation space +approxSpace = ApproximationSpace (dom) + +approxSpace:add_fct("u", "Lagrange",1,allSubsets) +approxSpace:add_fct("v", "Lagrange",1,allSubsets) +approxSpace:add_fct("p", "Lagrange",1,allSubsets) +approxSpace:add_fct("c", "Lagrange",1,allSubsets) + + +approxSpace:init_levels() +approxSpace:init_top_surface() +approxSpace:print_statistic() + +util.solver.defaults.approxSpace = approxSpace + + +------------------------------------------------------------------------------------------ +-- Compose the discretization +------------------------------------------------------------------------------------------ + + +-- Order the DoFs: +OrderLex (approxSpace, "x") +--OrderCuthillMcKee(approxSpace,true) +-- grid function for the solution +u = GridFunction(approxSpace) +u:set(0) + +-------------------------------- +-------------------------------- +-- Lua Functions +-------------------------------- +-------------------------------- + + +---------------------------------------------------------------------- Initial Velocity + +function ValueX(x,y) + hh=15 + return inflow*( (hh - y) * (4*y ) / (hh * hh)) +end + + +function ValueY(x,y) + return 0.0 +end +function ValueZ(x,y) + return 0 +end +---------------------------------------------------------------------- Initial Pressure + +function ValueP(x,y) + return 0.0 +end +-------------------------------------------------------------- Initial VolumeFraction + +function ValueC(x,y) + return 0.0 +end + +---------------------------------------------------------------------- Boundary Condition + +----------------------------------------------------------- Inlet +function BoundaryVolumeFraction(x,y) + return c_init +end + + +function inflowVel2d(x, y, t) + + return ValueX(x,y),ValueY(x,y) +end + + +------------------------------------------------------------------------------------------ +-- Compose the discretization of NavierStokes Equation +------------------------------------------------------------------------------------------ + +-- inner space +NavierStokesDisc = NavierStokesFV1(fct_cmp_tbl, {"Inner"}) + +NavierStokesDisc:set_peclet_blend(bPecletBlend) +NavierStokesDisc:set_exact_jacobian(bExactJac) +NavierStokesDisc:set_laplace(bNoLaplace) +NavierStokesDisc:set_stokes(bStokes) +NavierStokesDisc:set_upwind (upwind) +NavierStokesDisc:set_stabilization (stab, diffLength) +NavierStokesDisc:set_pac_upwind (bPac) + + +----------------------------------------------------------------------- Set inputs + +NavierStokesDisc:set_density(rho_a) +NavierStokesDisc:set_kinematic_viscosity(nu_a) + + + +------------------------------------------------- Set NS boundary conditions + +-- boundary condition at the inlet +InletDisc = NavierStokesInflow (NavierStokesDisc) +InletDisc:add ("inflowVel2d", "Left") + +-- boundary condition at the outlet +OutletDisc = NavierStokesNoNormalStressOutflow (NavierStokesDisc) +OutletDisc:add ("Right") + +-- boundary condition at the walls +WallDisc = NavierStokesWall (NavierStokesDisc) +WallDisc:add (" Bottom,Top") + +print("Navier Stokes Equation created.") + +------------------------------------------------------------------------------------------ +-- Compose the discretization of Transport Equation +------------------------------------------------------------------------------------------ + +TransportEq = ConvectionDiffusion("c", "Inner", "fv1") +TransportEq:set_velocity(NavierStokesDisc:velocity_ip()) +TransportEq:set_upwind(UpwindFV1(upwind)) + + +------------------------------------------------- Set boundary conditions + + +-- create dirichlet boundary for concentration +dirichletBND = DirichletBoundary() +dirichletBND:add(BoundaryVolumeFraction, "c", "Left") +dirichletBND:add(0.0, "c", "Top,Bottom") + + +OutflowBND = ConvectionDiffusionOutflowFV1(TransportEq) +OutflowBND:add( "Right") + + +print("Transport Equation created.") + + +--------------------------------------------------------------------------------------- +-- Domain Disc +--------------------------------------------------------------------------------------- + +domainDisc = DomainDiscretization(approxSpace) +domainDisc:add(NavierStokesDisc) +domainDisc:add(WallDisc) +domainDisc:add(InletDisc) +domainDisc:add(OutletDisc) + + +domainDisc:add(TransportEq) +domainDisc:add(OutflowBND) +domainDisc:add(dirichletBND) + +-- create operator from discretization + + +-- create time discretization +if timeMethod=="cn" then + timeDisc = ThetaTimeStep(domainDisc) + timeDisc:set_theta(0.5) -- Crank-Nicolson method +end +if timeMethod=="euler" then + timeDisc = ThetaTimeStep(domainDisc) + timeDisc:set_theta(1) -- implicit Euler +end +if timeMethod=="fracstep" then + timeDisc = ThetaTimeStep(domainDisc,"FracStep") +end +if timeMethod=="alex" then + timeDisc = ThetaTimeStep(domainDisc, "Alexander") +end + +-- create the assembled operator for the solver +op = AssembledOperator(timeDisc) + + +op:init() + +------------------------------------------------------------------------------------------ +-- Set up the solver +------------------------------------------------------------------------------------------ + +solverDesc = +{ + type = "newton", + linSolver = + { + type = "bicgstab", + precond = + { + type = "gmg", + rap = false, + smoother = + { + type = "ilu", + beta = ilu_beta, + consistentInterfaces = true + }, + preSmooth = 1, + postSmooth = 1, + baseSolver = "lu", + baseLevel = numPreRefs + }, + convCheck = + { + type = "standard", + iterations = 128, + absolute = 1e-12, + reduction = 0.5e-8, + verbose = true + } + }, + lineSearch = + { + type = "standard", + maxSteps = 3, + lambdaStart = 1, + lambdaReduce = 0.5, + acceptBest = true, + checkAll = false + }, + convCheck = + { + type = "standard", + iterations = 256, + absolute = 2.5e-12, + reduction = 1e-8, + verbose = true + } +} + +solver = util.solver.CreateSolver(solverDesc) + +------------------------------------------------------------------------------------------ +-- Prepare the initial values +------------------------------------------------------------------------------------------ + +-- grid function for the solution +u = GridFunction (approxSpace) + +-- set the initial guess +u:set (0) +Interpolate(0.0, u, "u") +Interpolate("ValueY", u, "v") +Interpolate("ValueP", u, "p") +Interpolate("ValueC", u, "c") + +-- order the DoFs: +--OrderCuthillMcKee (approxSpace, true) +OrderLex (approxSpace, "x") + + +------------------------------------------------------------------------------------------ +-- Apply the solver +------------------------------------------------------------------------------------------ +-- start +time = 0 +step = 0 + + +vtk_file_name = file_name .. "-lev" .. numRefs +if bStokes then + vtk_file_name = vtk_file_name .. "-Stokes" +end + -- write start solution + print("Writing start values") + out = VTKOutput() + out:clear_selection() + out:select_all(false) + out:select_nodal ("u,v", "velocity") + out:select_nodal ("u", "u") + out:select_nodal ("v", "v") + out:select_nodal ("p", "p") + out:select_nodal ("c", "c") + + + out:print_subsets(vtk_file_name, u,allSubsets,0,0) + +solver:init(op) + + +if solver:prepare(u) == false then + print ("Newton solver prepare failed.") exit() +end + +-- create new grid function for old value +uOld = u:clone() + +tBefore = os.clock() + +-- store grid function in vector of old solutions +solTimeSeries = SolutionTimeSeries() +solTimeSeries:push(uOld, time) + + +for step = 1, numTimeSteps do + print("++++++ TIMESTEP " .. step .. " BEGIN ++++++") + + -- choose time step + do_dt = dt + + -- setup time Disc for old solutions and timestep + timeDisc:prepare_step(solTimeSeries, do_dt) + + -- prepare newton solver + if solver:prepare(u) == false then + print ("Newton solver failed at step "..step.."."); exit(); + end + + if solver:apply(u) == false then + print ("Newton solver failed at step "..step.."."); exit(); + + end + + -- update new time + time= timeDisc:future_time() + + -- get oldest solution + oldestSol = solTimeSeries:oldest() + + -- copy values into oldest solution (we reuse the memory here) + VecAssign(oldestSol, u) + + -- push oldest solutions with new values to front, oldest sol pointer is poped from end + solTimeSeries:push_discard_oldest(oldestSol, time) + + + -- compute CFL number + + CFL=cflNumber(u,do_dt) + print("DT=" .. dt .. ""); + print("Time=" .. time .. ""); + + + if step % outputFactor == 0 then + + out = VTKOutput() + out:clear_selection() + out:select_all(false) + out:select_nodal("u,v", "velocity") + out:select_nodal("u", "u") + out:select_nodal("v", "v") + out:select_nodal("p", "p") + out:select_nodal("c", "c") + + + out:print_subsets(vtk_file_name, u,allSubsets,step,time) + print(" ") + end + print("++++++ TIMESTEP " .. step .. " END ++++++") +end + +tAfter = os.clock() +solver:print_average_convergence() +print("Computation took " .. tAfter-tBefore .. " seconds.") + +print("done.")