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

improved setup for simulation (ii) and (iii)

Browse source
This commit is contained in:
jakob.schratter 2026-01-27 09:14:25 +01:00
commit d8fc085d57
4 changed files with 5443 additions and 85 deletions
Download patch file Download diff file Expand all files Collapse all files

6
mgrid_2/getmatrix.cpp
View file

@ -637,8 +637,8 @@ void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
int const EdgeNode1 = RobinEdgeNodes[2*i]; int const EdgeNode1 = RobinEdgeNodes[2*i];
int const EdgeNode2 = RobinEdgeNodes[2*i + 1]; int const EdgeNode2 = RobinEdgeNodes[2*i + 1];
cout << "Edge number " << EdgeNumber << ", subdomain: " << subdomain << " "; // cout << "Edge number " << EdgeNumber << ", subdomain: " << subdomain << " ";
cout << "Node1: " << EdgeNode1 << " Node2: " << EdgeNode2 << " " << endl; // cout << "Node1: " << EdgeNode1 << " Node2: " << EdgeNode2 << " " << endl;
double x_1 = Coordinates[2*EdgeNode1]; double x_1 = Coordinates[2*EdgeNode1];
double y_1 = Coordinates[2*EdgeNode1 + 1]; double y_1 = Coordinates[2*EdgeNode1 + 1];
@ -653,7 +653,7 @@ void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
int ij = fetch(EdgeNode1, EdgeNode2); int ij = fetch(EdgeNode1, EdgeNode2);
int ji = fetch(EdgeNode2, EdgeNode1); int ji = fetch(EdgeNode2, EdgeNode1);
cout << "ii: " << ii << ", jj: " << jj << ", ij: " << ij << ", ji: " << ji << endl; //cout << "ii: " << ii << ", jj: " << jj << ", ij: " << ij << ", ji: " << ji << endl;
_sk[ii] += EdgeLength*alpha/3; _sk[ii] += EdgeLength*alpha/3;
_sk[jj] += EdgeLength*alpha/3; _sk[jj] += EdgeLength*alpha/3;

82
mgrid_2/main.cpp
View file

@ -7,6 +7,7 @@
#include <cassert> #include <cassert>
#include <chrono> // timing #include <chrono> // timing
#include <cmath> #include <cmath>
#include <fstream>
#include <iostream> #include <iostream>
#include <omp.h> #include <omp.h>
using namespace std; using namespace std;
@ -88,8 +89,8 @@ int main(int argc, char **argv )
double goal_perc = 60.0; double goal_perc = 60.0;
double time_count = 0; double time_count = 0;
while (average_cup_temperature < goal_temp) //while (average_cup_temperature < goal_temp)
// while (percentage_temp_reached < goal_perc) while (percentage_temp_reached < goal_perc)
//for (int step = 0; step < steps; ++step) //for (int step = 0; step < steps; ++step)
{ {
vector<double> G(Mdt.Nrows(), 0.0); vector<double> G(Mdt.Nrows(), 0.0);
@ -119,39 +120,58 @@ int main(int argc, char **argv )
mesh_c.Visualize(uv); mesh_c.Visualize(uv);
// // ################################## SIMULATION (ii) ##################################
// double u0_coffee = 85.0;
// mesh_c.Init_Solution_mult(uv, 1, [u0_coffee](double x, double y) -> double { return u0_coffee; }); // fluid
// t3 = system_clock::now(); // start timer // output vector for simulation (iii)
// double average_coffee_temperature = u0_coffee; ofstream output_file("../solid-cpp/uv_1.txt");
for (double node_value : uv)
// time_count = 0; {
// while (average_coffee_temperature > 50.0) output_file << node_value << endl;
// { }
// vector<double> G(Mdt.Nrows(), 0.0); output_file.close();
// Mdt.Mult(G, uv); // G = M/dt * u_{n}
// vector<double> H = fv;
// for (size_t i = 0; i < Mdt.Nrows(); ++i)
// {
// H[i] += G[i]; // H = F + G
// }
// JacobiSolve(SK, H, uv); // solve: (M/dt + K + C) * u_{n+1} = F + M/dt * u_{n}
// // ----- SK ----- ------ H -------
// average_coffee_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 1);
// cout << "Average coffee temperature: " << average_coffee_temperature << " after " << time_count << " seconds. " << endl;
// time_count += dt;
// }
// t4 = system_clock::now(); // stop timer
// duration = duration_cast<microseconds>(t4 - t3); // duration in microseconds // ################################## SIMULATION (ii) ##################################
// t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds double u0_coffee = 85.0;
// cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl; mesh_c.Init_Solution_mult(uv, 1, [u0_coffee](double x, double y) -> double { return u0_coffee; }); // fluid
mesh_c.Visualize(uv);
t3 = system_clock::now(); // start timer
double average_coffee_temperature = u0_coffee;
goal_temp = 50.0;
goal_perc = 60.0;
percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 1, goal_temp);
time_count = 0;
//while (average_coffee_temperature > goal_temp)
while (percentage_temp_reached > goal_perc)
{
vector<double> G(Mdt.Nrows(), 0.0);
Mdt.Mult(G, uv); // G = M/dt * u_{n}
vector<double> H = fv;
for (int i = 0; i < Mdt.Nrows(); ++i)
{
H[i] += G[i]; // H = F + G
}
JacobiSolve(SK, H, uv); // solve: (M/dt + K + C) * u_{n+1} = F + M/dt * u_{n}
// ----- SK ----- ------ H -------
average_coffee_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 1);
percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 1, goal_temp);
cout << "Average coffee temperature: " << average_coffee_temperature << " after " << time_count << " seconds. " << endl;
cout << "% of elements above temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl;
time_count += dt;
}
t4 = system_clock::now(); // stop timer
duration = duration_cast<microseconds>(t4 - t3); // duration in microseconds
t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl;
mesh_c.Visualize(uv);
return 0; return 0;

83
solid-cpp/ownSolver.cpp
View file

@ -7,6 +7,7 @@
#include <cassert> #include <cassert>
#include <chrono> // timing #include <chrono> // timing
#include <cmath> #include <cmath>
#include <fstream>
#include <iostream> #include <iostream>
#include <omp.h> #include <omp.h>
@ -76,61 +77,26 @@ int main(int argc, char **argv )
// ########################################## // ##########################################
// ################################## SIMULATION (i) ##################################
// Initialize temperature u_0
vector<double> uv(SK.Nrows(), 0.0); // temperature
mesh_c.Init_Solution_mult(uv, 0, [u0_mug](double x, double y) -> double { return u0_mug; }); // mug
mesh_c.Init_Solution_mult(uv, 1, [u0_fluid](double x, double y) -> double { return u0_fluid; }); // fluid
mesh_c.Init_Solution_mult(uv, 2, [u0_air](double x, double y) -> double { return u0_air; }); // air
//mesh_c.Visualize(uv);
auto t3 = system_clock::now(); // start timer
double average_cup_temperature = u0_mug;
double percentage_temp_reached = 0.0;
double goal_temp = 67.0;
double goal_perc = 60.0;
double time_count = 0;
while (average_cup_temperature < goal_temp)
// while (percentage_temp_reached < goal_perc)
//for (int step = 0; step < steps; ++step)
{
vector<double> G(Mdt.Nrows(), 0.0);
Mdt.Mult(G, uv); // G = M/dt * u_{n}
vector<double> H = fv;
for (size_t i = 0; i < Mdt.Nrows(); ++i)
{
H[i] += G[i]; // H = F + G
}
JacobiSolve(SK, H, uv); // solve: (M/dt + K + C) * u_{n+1} = F + M/dt * u_{n}
// ----- SK ----- ------ H -------
average_cup_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 0);
percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 0, goal_temp);
cout << "Average cup temperature: " << average_cup_temperature << " after " << time_count << " seconds. " << endl;
cout << "% of elements reached temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl;
time_count += dt;
}
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 << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl;
mesh_c.Visualize(uv);
// ################################## SIMULATION (iii) ################################## // ################################## SIMULATION (iii) ##################################
// read vector from simulation (i)
vector<double> uv(SK.Nrows(), 0.0);
ifstream input_file("uv_1.txt");
for (size_t i = 0; i < uv.size(); ++i)
{
input_file >> uv[i];
}
double u0_coffee = 85.0; double u0_coffee = 85.0;
mesh_c.Init_Solution_mult(uv, 1, [u0_coffee](double x, double y) -> double { return u0_coffee; }); // fluid mesh_c.Init_Solution_mult(uv, 1, [u0_coffee](double x, double y) -> double { return u0_coffee; }); // fluid
t3 = system_clock::now(); // start timer mesh_c.Visualize(uv);
auto t3 = system_clock::now(); // start timer
double average_coffee_temperature = u0_coffee; double average_coffee_temperature = u0_coffee;
percentage_temp_reached = 0.0; double goal_temp = 50.0;
double goal_perc = 60.0;
double percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 1, goal_temp);
// ------------------------ initialize preCICE ------------------------ // ------------------------ initialize preCICE ------------------------
int commRank = 0; int commRank = 0;
@ -172,9 +138,9 @@ int main(int argc, char **argv )
goal_temp = 50.0; goal_temp = 50.0;
goal_perc = 60.0; goal_perc = 60.0;
time_count = 0; double time_count = 0;
while (average_cup_temperature > goal_temp) //while (average_coffee_temperature > goal_temp)
// while (percentage_temp_reached > goal_perc) while (percentage_temp_reached > goal_perc)
{ {
preciceDt = participantSolid.getMaxTimeStepSize(); preciceDt = participantSolid.getMaxTimeStepSize();
solverDt = 1.0; solverDt = 1.0;
@ -201,9 +167,10 @@ int main(int argc, char **argv )
// ----- SK ----- ------ H ------- // ----- SK ----- ------ H -------
average_coffee_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 1); average_coffee_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 1);
percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 0, goal_temp); percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 1, goal_temp);
cout << "Average coffee temperature: " << average_coffee_temperature << " after " << time_count << " seconds. " << endl; cout << "Average coffee temperature: " << average_coffee_temperature << " after " << time_count << " seconds. " << endl;
cout << "% of elements reached temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl; cout << "% of elements above temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl;
time_count += dt;
// ----- write the heat-flux, so openFOAM can read it // ----- write the heat-flux, so openFOAM can read it
{ {
@ -218,9 +185,9 @@ int main(int argc, char **argv )
time_count += dt; time_count += dt;
} }
t4 = system_clock::now(); // stop timer auto t4 = system_clock::now(); // stop timer
duration = duration_cast<microseconds>(t4 - t3); // duration in microseconds auto duration = duration_cast<microseconds>(t4 - t3); // duration in microseconds
t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl; cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl;

5371
solid-cpp/uv_1.txt Normal file
View file

File diff suppressed because it is too large Load diff