diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp index d20e630..ab55f9a 100644 --- a/mgrid_2/getmatrix.cpp +++ b/mgrid_2/getmatrix.cpp @@ -383,7 +383,7 @@ void FEM_Matrix::Derive_Matrix_Pattern_slow() return; } -void FEM_Matrix::CalculateLaplaceMult(vector &f) +void FEM_Matrix::CalculateLaplace_mult(vector &f) { cout << "\n############ FEM_Matrix::CalculateLaplaceMult "; double tstart = omp_get_wtime(); // OpenMP @@ -394,9 +394,6 @@ void FEM_Matrix::CalculateLaplaceMult(vector &f) for (int k = 0; k < _nrows; ++k) { _sk[k] = 0.0; } - for (int k = 0; k < _nrows; ++k) { - f[k] = 0.0; - } double ske[3][3], fe[3]; // Loop over all elements diff --git a/mgrid_2/getmatrix.h b/mgrid_2/getmatrix.h index 39cef6f..d3bd93d 100644 --- a/mgrid_2/getmatrix.h +++ b/mgrid_2/getmatrix.h @@ -410,8 +410,8 @@ class FEM_Matrix: public CRS_Matrix * The penalty method * is used for incorporating the given values @p u. * - * @param[in] u_out outside temperature * @param[in,out] f load vector + * @param[in] u_out outside temperature */ void ApplyRobinBC_mult(std::vector &f, const double u_out); @@ -713,6 +713,7 @@ void CalcElem_Masse(int const ial[3], double const xc[], double ske[3][3]); * for specific volumetric heat capacity in subdomain. * @param[in] ial node indices of the three element vertices * @param[in] xc vector of node coordinates with x(2*k,2*k+1) as coordinates of node k + * @param[in] c volumetric heat capacity of element * @param[out] ske element stiffness matrix */ void CalcElem_MasseSpecific(int const ial[3], double const xc[], double const c, double ske[3][3]); diff --git a/mgrid_2/main.cpp b/mgrid_2/main.cpp index 5414edc..7dddf2a 100644 --- a/mgrid_2/main.cpp +++ b/mgrid_2/main.cpp @@ -13,17 +13,8 @@ using namespace std; using namespace std::chrono; // timing -void Test_solver(int lev, vector const &nthreads, Multigrid &ggm); - int main(int argc, char **argv ) { -#undef MG -#ifndef MG - // Jacobi iteration - //int nrefine = 0; - //if (argc > 1) nrefine = atoi(argv[1]); - - // generating the mesh Mesh const mesh_c("../generate_mesh/coffee_cup.txt", "../generate_mesh/coffee_cup_sd.txt"); bool ba = mesh_c.checkObtuseAngles(); @@ -35,28 +26,31 @@ int main(int argc, char **argv ) //mesh.Debug(); //mesh_c.DebugEdgeBased(); + // ########################################## + // Assembling + // ########################################## // Initializing FEM matrix !pattern! (only zero entries now) - FEM_Matrix SK(mesh_c); // CRS matrix + FEM_Matrix SK(mesh_c); // CRS matrix //SK.writeBinary("sparseMatrix.bin"); //SK.Debug(); + vector fv(SK.Nrows(), 0.0); + SK.CalculateRHS(fv, [](double x, double y) {return 0;}); // r.h.s. - // Initialize RHS - vector fv(SK.Nrows(), 0.0); // r.h.s. + SK.CalculateLaplace_mult(fv); // stiffness matrix + SK.AddMass_mult(fv); // mass matrix + SK.ApplyRobinBC_mult(fv, 18.0); // apply Robin bnd - // Calculate stiffness matrix entries - SK.CalculateLaplaceMult(fv); // matrix - //SK.Debug(); + SK.CheckRowSum(); + SK.CheckMatrix(); - // Add mass matrix entries - SK.AddMass_mult(fv); - - // Calculate RHS - SK.CalculateRHS(fv, [](double x, double y) {return 0;}); - //SK.CheckRowSum(); - //SK.CheckMatrix(); + // ########################################## + // Timestepping + // ########################################## + double dt = 1.0; // time step + int steps = 100; // number of time iterations // Initialize temperature vector uv(SK.Nrows(), 0.0); // temperature @@ -64,88 +58,21 @@ int main(int argc, char **argv ) mesh_c.Init_Solution_mult(uv, 1, [](double x, double y) -> double { return 80; }); // fluid mesh_c.Init_Solution_mult(uv, 2, [](double x, double y) -> double { return 18; }); // air - // Apply BC - SK.ApplyRobinBC_mult(fv, 18.0); - - - // Solve - auto exact_sol(uv); - //SK.Mult(fv,uv); + for (int i = 0; i < steps; ++i) + { + // TODO + } auto t3 = system_clock::now(); // start timer - //JacobiSolve(SK, fv, uv ); // solve the system of equations - auto t4 = system_clock::now(); // stop timer + auto duration = duration_cast(t4 - t3); // duration in microseconds double t_diff = static_cast(duration.count()) / 1e6; // overall duration in seconds cout << "JacobiSolve: timing in sec. : " << t_diff << endl; - // Calculate error and visualize - auto [val, idx] = findLargestAbsError(exact_sol, uv, 1e+6, 100); - - //mesh.Visualize(getAbsError(exact_sol, uv)); mesh_c.Visualize(uv); -#else - // multigrid iteration - int nrefine = 3; - if (argc > 1) nrefine = atoi(argv[1]); - - //Multigrid ggm(Mesh("square_tiny.txt"),nrefine); - Multigrid ggm(Mesh("square_100.txt"), nrefine); - - ggm.DefineOperators(); - - cout << "\n############# SOLVE ###############\n"; - - double tstart = omp_get_wtime(); // OpenMP - - //ggm.JacobiSolve(my_level); - //ggm.MG_Step(my_level, 1, true, 1); - ggm.MG_Solve(2, 1e-6, 1); - - double t1 = omp_get_wtime() - tstart; // OpenMP - cout << "MgSolve: timing in sec. : " << t1 << " for " << ggm.Ndofs() << " dofs" << endl; - - - Test_solver(nrefine - 1, {1, 2, 4, 8, 16, 32, 64, 128, 256}, ggm); - - //int my_level=nrefine-1; - //const auto &ml=ggm.GetMesh(my_level); - //const auto &sl=ggm.GetSolution(my_level); - //ml.Visualize(sl); - //////ml.Visualize_paraview(sl); - ////ml.Export_scicomp("level_"+to_string(my_level)); - - //int my_level=nrefine-1; - //const auto &mesh=ggm.GetMesh(my_level); - //const auto &uv=ggm.GetSolution(my_level); - //vector exact_sol(size(uv)); - //mesh.SetValues(exact_sol, [](double x, double y) -> double { - //return x * x * std::sin(2.5 * M_PI * y); - //} ); - //mesh.Visualize(getAbsError(exact_sol, uv)); - -#endif return 0; -} - - -void Test_solver(int /*lev*/, vector const &nthreads, Multigrid &ggm) -{ - cout << endl << endl << "-------------------------------------" << endl; - cout << "MgSolve: timing in sec. for " << ggm.Ndofs() << " dofs" << endl; - cout << "sec threads" << endl; - vector mg_time(size(nthreads), -1.0); - - for (size_t k = 0; k < size(nthreads); ++k) { - omp_set_num_threads(nthreads.at(k)); - double tstart = omp_get_wtime(); - ggm.MG_Solve(2, 1e-6, 1); - double t1 = omp_get_wtime() - tstart; - mg_time.at(k) = t1; - cout << t1 << " : " << nthreads.at(k) << endl; - } -} +} \ No newline at end of file