From c9deacfd3c449d3abe0dcc662a7a24273bfe4a3e Mon Sep 17 00:00:00 2001 From: "jakob.schratter" Date: Mon, 26 Jan 2026 16:16:12 +0100 Subject: [PATCH] C++ solver part --- solid-cpp/CMakeLists.txt | 15 +++ solid-cpp/clean.sh | 7 ++ solid-cpp/ownSolver.cpp | 220 +++++++++++++++++++++++++++++++++++++++ solid-cpp/run.sh | 14 +++ 4 files changed, 256 insertions(+) create mode 100644 solid-cpp/CMakeLists.txt create mode 100755 solid-cpp/clean.sh create mode 100644 solid-cpp/ownSolver.cpp create mode 100755 solid-cpp/run.sh diff --git a/solid-cpp/CMakeLists.txt b/solid-cpp/CMakeLists.txt new file mode 100644 index 0000000..9b025a6 --- /dev/null +++ b/solid-cpp/CMakeLists.txt @@ -0,0 +1,15 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 3.10.2) +SET(TARGET "ownSolver") +PROJECT(${TARGET} LANGUAGES CXX DESCRIPTION "ownSolver") + +# The ownSolver requires c++14 +SET(CMAKE_CXX_STANDARD 14) +SET(CMAKE_CXX_EXTENSIONS OFF) +SET(CMAKE_CXX_STANDARD_REQUIRED ON) + +FIND_PACKAGE(precice 3.0 REQUIRED CONFIG) +ADD_EXECUTABLE( + ${TARGET} + ${TARGET}.cpp) + +TARGET_LINK_LIBRARIES(${TARGET} PRIVATE precice::precice) diff --git a/solid-cpp/clean.sh b/solid-cpp/clean.sh new file mode 100755 index 0000000..e8bdcf6 --- /dev/null +++ b/solid-cpp/clean.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +. ../preCICE_tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/solid-cpp/ownSolver.cpp b/solid-cpp/ownSolver.cpp new file mode 100644 index 0000000..6ecfe78 --- /dev/null +++ b/solid-cpp/ownSolver.cpp @@ -0,0 +1,220 @@ +#include "geom.h" +#include "getmatrix.h" +#include "jacsolve.h" +#include "userset.h" +#include "vdop.h" + +#include +#include // timing +#include +#include +#include + +#include "precice/precice.hpp" +#include +using namespace precice; + +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 = 100; // 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 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} + // ########################################## + + + // ################################## SIMULATION (i) ################################## + + // Initialize temperature u_0 + vector 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 time_count = 0; + while (average_cup_temperature < 67.0) + //for (int step = 0; step < steps; ++step) + { + vector G(Mdt.Nrows(), 0.0); + Mdt.Mult(G, uv); // G = M/dt * u_{n} + + vector 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 << " after " << time_count << " seconds. " << endl; + time_count += dt; + } + 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 << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl; + + mesh_c.Visualize(uv); + + // ################################## SIMULATION (iii) ################################## + 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 + double average_coffee_temperature = u0_coffee; + + // ------------------------ initialize preCICE ------------------------ + int commRank = 0; + int commSize = 1; + std::string configFileName("precise-config.xml"); + std::string solverName("Solid"); + std::string meshName("Solid-Mesh"); + + std::cout << "Running Solid-solver with preCICE config file \"" << configFileName << "\" and participant name \"" << solverName << "\".\n"; + + Participant participantSolid(solverName, configFileName, commRank, commSize); + + int meshDim = participantSolid.getMeshDimensions(meshName); // gets mesh dimensions (=2) from config file + int numberOfVertices; // number of vertices at "wet" surface + { + // TODO: DETERMINE NUMBER OF VERTICES AT TOP + } + std::vector coords(numberOfVertices*meshDim); // coords of vertices at "wet" surface + { + // TODO: DETERMINE COORDINATES OF VERTICES AT TOP + } + std::vector vertexIDs(numberOfVertices); + + participantSolid.setMeshVertices(meshName, coords, vertexIDs); + + // ----- initialize read- and write-data + std::vector temperature(numberOfVertices * meshDim); // read-data + std::vector heatFlux(numberOfVertices * meshDim); // write-data + + + double solverDt = 1.0; // solver time step size + double preciceDt = 1.0; // maximum precice time step size + dt = 1.0; // actual time step size + + participantSolid.initialize(); + + + // ------------------------ timestepping ------------------------ + time_count = 0; + while (average_coffee_temperature > 50.0) + { + preciceDt = participantSolid.getMaxTimeStepSize(); + solverDt = 1.0; + dt = min(preciceDt, solverDt); + + // ----- read temperature computed by openFOAM simulation ----- + participantSolid.readData("Solid-Mesh", "Temperature", vertexIDs, dt, temperature); + { + // TODO: SET THE MESH TO THE INCOMING TEMPERATURE + // "initNewTemp(temperature)" + } + + // ----- solve time step ----- + vector G(Mdt.Nrows(), 0.0); + Mdt.Mult(G, uv); // G = M/dt * u_{n} + + vector 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; + + // ----- write the heat-flux, so openFOAM can read it + { + // TODO: COMPUTE THE HEAT-FLUX FROM THE NEW uv TEMPERATURE SOLUTION + // "computeHeatFlux(heatFlux)" + } + participantSolid.writeData("Solid-Mesh", "Heat-Flux", vertexIDs, heatFlux); + + // ----- advance preCICE ----- + participantSolid.advance(dt); + + + time_count += dt; + } + t4 = system_clock::now(); // stop timer + duration = duration_cast(t4 - t3); // duration in microseconds + t_diff = static_cast(duration.count()) / 1e6; // overall duration in seconds + cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl; + + + // ----- free data, close communication ----- + participantSolid.finalize(); + + + + return 0; +} \ No newline at end of file diff --git a/solid-cpp/run.sh b/solid-cpp/run.sh new file mode 100755 index 0000000..ec2177d --- /dev/null +++ b/solid-cpp/run.sh @@ -0,0 +1,14 @@ +#!/bin/bash +set -e -u + +. ../preCICE_tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +solver=./ownSolver +if [ -f "${solver}" ]; then + ${solver} +else + echo "Unable to locate the executable ${solver}. Have a look at the README for building instructions." +fi + +close_log