/*********************************************************************************** OpenCS Project: www.daetools.com Copyright (C) Dragan Nikolic ************************************************************************************ OpenCS is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License version 3 as published by the Free Software Foundation. OpenCS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the OpenCS software; if not, see <http://www.gnu.org/licenses/>. ***********************************************************************************/ #include <mpi.h> #include "burgers_equations_cv_2d.h" #include <OpenCS/models/cs_model_builder.h> #include <OpenCS/models/cs_partitioners.h> #include <OpenCS/simulators/cs_simulators.h> using namespace cs; /* The model describes a viscous fluid flow defined by the Burgers equations * (transient convection-diffusion equations): * du/dt + (d(u.u)/dx + d(u.v)/dy) - ni * (d2u/dx2 + d2u/dy2) = Su(x,y) * dv/dt + (d(v.u)/dx + d(v.v)/dy) - ni * (d2v/dx2 + d2v/dy2) = Sv(x,y) * * presented in Appendix C of: * - K. Salari, P. Knupp. Code Verification by the Method of Manufactured Solutions. * SAND2000 – 1444 (2000), <https://doi.org/10.2172/759450>. * and represents of re-implementation of DAE Tools code verifications tutorial CV-4: * docs/tutorials-cv.html#tutorial-cv-4. * This problem uses a formal code verification technique: * the Method of Manufactured Solutions for checking the correctness of the numerical software. * The procedure is the following: * (1) a function representing the solution is selected (the manufactured solution), * (2) the new source term for the original problem is computed, and * (3) the original problem is solved with the new source term. * The computed numerical solution should be equal to the manufactured one. * The manufactured solutions are given by: * um(x,y) = u0 * (sin(x^2 + y^2 + w0*t) + eps) * vm(x,y) = v0 * (cos(x^2 + y^2 + w0*t) + eps) * * The terms in the new sources Su and Sv are computed by differentiation of * the manufactured solution equations for um and vm. * Su(x,y) = dum/dt + (d(um.um)/dx + d(um.vm)/dy) - ni * (d2um/dx2 + d2um/dy2) * Sv(x,y) = dvm/dt + (d(vm.um)/dx + d(vm.vm)/dy) - ni * (d2vm/dx2 + d2vm/dy2) * * The Dirichlet boundary conditions are applied at boundaries: * u(x, y) = um(x, y) * v(x, y) = vm(x, y) * for x = {0, Nx-1} and y = {0, Ny-1}. * * The equations are distributed on a non-uniform structured rectangular 2D grid * (-0.1, 0.7)x(0.2, 0.8) with 10x8, 20x16, 40x32 and 80x64 grid points to avoid symmetry * for a more thorough check of the computed solution. */ int main(int argc, char *argv[]) { int Nexps = 4; int Nxs[] = {10, 20, 40, 80}; int Nys[] = {8, 16, 32, 64}; MPI_Init(&argc, &argv); for(int i = 0; i < Nexps; i++) { uint32_t Nx = Nxs[i]; uint32_t Ny = Nys[i]; printf("########################################################\n"); printf(" Burgers equations (groups) in 2D: %u x %u grid\n", Nx, Ny); printf("########################################################\n"); burgers_groups_2D burgers(Nx, Ny); uint32_t Ndofs = 0; uint32_t Nvariables = burgers.Nequations; real_t defaultVariableValue = 0.0; real_t defaultVariableTimeDerivative = 0.0; real_t defaultAbsoluteTolerance = 1e-6; std::string defaultVariableName = "x"; /* 1. Initialise model builder with the number of variables/equations. */ csModelBuilder_t mb; mb.Initialize_DAE_System(Nvariables, Ndofs, defaultVariableValue, defaultVariableTimeDerivative, defaultAbsoluteTolerance, defaultVariableName); printf("Model builder initialised\n"); /* 2. Create and set model equations using the provided time/variable/derivative/dof objects. */ const csNumber_t& time = mb.GetTime(); const std::vector<csNumber_t>& x_vars = mb.GetVariables(); const std::vector<csNumber_t>& xt_vars = mb.GetTimeDerivatives(); const std::vector<csNumber_t>& y_vars = mb.GetDegreesOfFreedom(); std::vector<csKernelPtr> kernels; std::vector<csEquation_t> equations; burgers.CreateEquations(time, x_vars, xt_vars, equations); printf("Model equations generated\n"); mb.SetModelEquations(equations, kernels); printf("Model equations set\n"); // Set variable names std::vector<std::string> names; burgers.GetVariableNames(names); mb.SetVariableNames(names); /* 3. Export model(s) to a specified directory for sequential and parallel simulations. */ // Set initial conditions. std::vector<real_t> x_0(Nvariables, 0.0); burgers.SetInitialConditions(x_0); mb.SetVariableValues(x_0); // Set the simulation options. csSimulationOptionsPtr options = mb.GetSimulationOptions(); options->LoadString(simulation_options); options->SetDouble("Simulation.TimeHorizon", 90.0); options->SetDouble("Simulation.ReportingInterval", 2.0); options->SetDouble("Solver.Parameters.RelativeTolerance", 1e-5); std::string outDir = std::to_string(Nx) + "x" + std::to_string(Ny); options->SetString("Simulation.OutputDirectory", outDir); std::vector<std::string> grps_0 = {"BoundaryConditions"}; std::vector<std::string> grps_1 = {"Burgers_u", "Burgers_v", "Burgers_um", "Burgers_vm"}; /* (1) Sequential evaluator for boundary conditions. */ //options->SetString("Model.Evaluators.Device_0.Library", "Sequential"); //options->SetString("Model.Evaluators.Device_0.Name", "seq"); //options->SetStringArray("Model.Evaluators.Device_0.Groups", grps_0); /* (2a) OpenMP evaluator for u, v, um, vm groups. */ //options->SetString("Model.Evaluators.Device_1.Library", "OpenMP"); //options->SetString("Model.Evaluators.Device_1.Name", "openmp"); //options->SetInteger("Model.Evaluators.Device_1.Parameters.numThreads", 0); //options->SetStringArray("Model.Evaluators.Device_1.Groups", grps_1); // For Ncpu = 1: k = 3, rho = 1.0, alpha = 1e-1, w = 0.5 //options->SetInteger("LinearSolver.Preconditioner.Parameters.fact: level-of-fill", 3); //options->SetDouble ("LinearSolver.Preconditioner.Parameters.fact: relax value", 0.5); //options->SetDouble ("LinearSolver.Preconditioner.Parameters.fact: absolute threshold", 1e-1); //options->SetDouble ("LinearSolver.Preconditioner.Parameters.fact: relative threshold", 1.0); std::string simulationOptions_seq = options->ToString(); // Generate a single model (no graph partitioner required). std::string inputFilesDirectory = "dae_example_4_cv-sequential"; { uint32_t Npe = 1; std::vector<std::string> balancingConstraints; csGraphPartitioner_t* gp = NULL; std::vector<csModelPtr> models_sequential = mb.PartitionSystem(Npe, gp, balancingConstraints, true); csModelBuilder_t::ExportModels(models_sequential, inputFilesDirectory, simulationOptions_seq); printf("Generated model for Npe = 1\n"); } /* 4. Simulate the sequential model using the libOpenCS_Simulators. */ printf("Simulation of '%s' (using libOpenCS_Simulators)\n\n", inputFilesDirectory.c_str()); csSimulate(inputFilesDirectory); } MPI_Finalize(); std::system("python dae_example_4_cv_plots.py"); return 0; }