jakob.schratter/SciFEM-Project_CoffeeMugSimulation_Celebic-Schratter
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

cleanup

Browse source
This commit is contained in:
dino.celebic 2026-01-24 21:24:09 +01:00
commit c53335074d
3 changed files with 25 additions and 100 deletions
Download patch file Download diff file Expand all files Collapse all files

5
mgrid_2/getmatrix.cpp
View file

@ -383,7 +383,7 @@ void FEM_Matrix::Derive_Matrix_Pattern_slow()
return;
}
void FEM_Matrix::CalculateLaplaceMult(vector<double> &f)
void FEM_Matrix::CalculateLaplace_mult(vector<double> &f)
{
cout << "\n############ FEM_Matrix::CalculateLaplaceMult ";
double tstart = omp_get_wtime(); // OpenMP
@ -394,9 +394,6 @@ void FEM_Matrix::CalculateLaplaceMult(vector<double> &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

3
mgrid_2/getmatrix.h
View file

@ -410,8 +410,8 @@ class FEM_Matrix: public CRS_Matrix
* The <a href="https://www.jstor.org/stable/2005611?seq=1#metadata_info_tab_contents">penalty method</a>
* 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<double> &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]);

117
mgrid_2/main.cpp
View file

@ -13,17 +13,8 @@ using namespace std;
using namespace std::chrono; // timing
void Test_solver(int lev, vector<int> 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<double> fv(SK.Nrows(), 0.0);
SK.CalculateRHS(fv, [](double x, double y) {return 0;}); // r.h.s.
// Initialize RHS
vector<double> 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<double> 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<microseconds>(t4 - t3); // duration in microseconds
double t_diff = static_cast<double>(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<double> 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<int> const &nthreads, Multigrid &ggm)
{
cout << endl << endl << "-------------------------------------" << endl;
cout << "MgSolve: timing in sec. for " << ggm.Ndofs() << " dofs" << endl;
cout << "sec threads" << endl;
vector<double> 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;
}
}
}