// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. // 3D flow over a cylinder benchmark example #include "navier_solver.hpp" #include using namespace mfem; using namespace navier; struct s_NavierContext { int order = 4; double kin_vis = 0.001; double t_final = 8.0; double dt = 1e-3; } ctx; void vel(const Vector &x, double t, Vector &u) { double xi = x(0); double yi = x(1); double zi = x(2); double U = 2.25; if (xi <= 1e-8) { u(0) = 16.0 * U * yi * zi * sin(M_PI * t / 8.0) * (0.41 - yi) * (0.41 - zi) / pow(0.41, 4.0); } else { u(0) = 0.0; } u(1) = 0.0; u(2) = 0.0; } int main(int argc, char *argv[]) { Mpi::Init(argc, argv); Hypre::Init(); int serial_refinements = 0; Mesh *mesh = new Mesh("box-cylinder.mesh"); for (int i = 0; i < serial_refinements; ++i) { mesh->UniformRefinement(); } if (Mpi::Root()) { std::cout << "Number of elements: " << mesh->GetNE() << std::endl; } auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh); delete mesh; // Create the flow solver. NavierSolver flowsolver(pmesh, ctx.order, ctx.kin_vis); flowsolver.EnablePA(true); // Set the initial condition. ParGridFunction *u_ic = flowsolver.GetCurrentVelocity(); VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel); u_ic->ProjectCoefficient(u_excoeff); // Add Dirichlet boundary conditions to velocity space restricted to // selected attributes on the mesh. Array attr(pmesh->bdr_attributes.Max()); // Inlet is attribute 1. attr[0] = 1; // Walls is attribute 3. attr[2] = 1; flowsolver.AddVelDirichletBC(vel, attr); double t = 0.0; double dt = ctx.dt; double t_final = ctx.t_final; bool last_step = false; flowsolver.Setup(dt); ParGridFunction *u_gf = flowsolver.GetCurrentVelocity(); ParGridFunction *p_gf = flowsolver.GetCurrentPressure(); ParaViewDataCollection pvdc("3dfoc", pmesh); pvdc.SetDataFormat(VTKFormat::BINARY32); pvdc.SetHighOrderOutput(true); pvdc.SetLevelsOfDetail(ctx.order); pvdc.SetCycle(0); pvdc.SetTime(t); pvdc.RegisterField("velocity", u_gf); pvdc.RegisterField("pressure", p_gf); pvdc.Save(); for (int step = 0; !last_step; ++step) { if (t + dt >= t_final - dt / 2) { last_step = true; } flowsolver.Step(t, dt, step); if (step % 10 == 0) { pvdc.SetCycle(step); pvdc.SetTime(t); pvdc.Save(); } if (Mpi::Root()) { printf("%11s %11s\n", "Time", "dt"); printf("%.5E %.5E\n", t, dt); fflush(stdout); } } flowsolver.PrintTimingData(); delete pmesh; return 0; }