112 lines
No EOL
3.9 KiB
C++
112 lines
No EOL
3.9 KiB
C++
#include "geom.h"
|
|
#include "getmatrix.h"
|
|
#include "jacsolve.h"
|
|
#include "userset.h"
|
|
#include "vdop.h"
|
|
|
|
#include <cassert>
|
|
#include <chrono> // timing
|
|
#include <cmath>
|
|
#include <iostream>
|
|
#include <omp.h>
|
|
using namespace std;
|
|
using namespace std::chrono; // timing
|
|
|
|
|
|
int main(int argc, char **argv )
|
|
{
|
|
// generating the mesh
|
|
Mesh const mesh_c("../generate_mesh/coffee_cup.txt", "../generate_mesh/coffee_cup_sd.txt");
|
|
bool ba = mesh_c.checkObtuseAngles();
|
|
if (ba) cout << "mesh corrected" << endl;
|
|
//mesh_c.DebugEdgeBased();
|
|
|
|
//gMesh_Hierarchy ggm(mesh_c, nrefine);
|
|
//const Mesh &mesh = ggm.finest();
|
|
//mesh.Debug();
|
|
//mesh_c.DebugEdgeBased();
|
|
|
|
// ##########################################
|
|
// Parameteres
|
|
// ##########################################
|
|
|
|
double dt = 1.0; // time step
|
|
//int steps = 200; // number of time iterations
|
|
|
|
double u0_mug = 18.0;
|
|
double u0_fluid = 80.0;
|
|
double u0_air = 18.0;
|
|
double u_out = 18.0;
|
|
|
|
// ##########################################
|
|
// Assembling
|
|
// ##########################################
|
|
|
|
// Initializing FEM matrix !pattern! (only zero entries now)
|
|
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;}); // rhs (f=0)
|
|
|
|
SK.CalculateLaplace_mult(fv); // stiffness matrix (+K)
|
|
SK.AddMass_mult(fv, 1.0/dt); // add mass matrix (+M/dt)
|
|
SK.ApplyRobinBC_mult(fv, u_out); // apply Robin-bnd (+C = +F)
|
|
// SK = M/dt + K + C = F
|
|
|
|
// SK.Debug();
|
|
// SK.CheckRowSum();
|
|
// SK.CheckMatrix();
|
|
|
|
FEM_Matrix Mdt(mesh_c);
|
|
Mdt.AddMass_mult(fv, 1.0/dt); // mass matrix (Mdt = M/dt)
|
|
|
|
// Mdt.Debug();
|
|
// Mdt.CheckRowSum();
|
|
// Mdt.CheckMatrix();
|
|
|
|
// ##########################################
|
|
// Timestepping (M/dt + K + C) * u_{n+1} = F + M/dt * u_{n}
|
|
// ##########################################
|
|
|
|
// 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;
|
|
while (average_cup_temperature < 67)
|
|
//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);
|
|
cout << "Average cup temperature: " << average_cup_temperature << endl;
|
|
}
|
|
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);
|
|
|
|
return 0;
|
|
} |