C++ solver part

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)

solid-cpp/clean.sh

@@ -0,0 +1,7 @@
+#!/bin/sh
+set -e -u
+
+. ../preCICE_tools/cleaning-tools.sh
+
+clean_precice_logs .
+clean_case_logs .

solid-cpp/ownSolver.cpp

@@ -0,0 +1,220 @@
+#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>
+
+#include "precice/precice.hpp"
+#include <sstream>
+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<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}
+    // ##########################################
+
+
+    // ################################## 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 time_count = 0;
+    while (average_cup_temperature < 67.0)
+    //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 << " after " << time_count << " seconds. " << 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) ##################################
+    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<double> coords(numberOfVertices*meshDim); // coords of vertices at "wet" surface
+    {
+        // TODO: DETERMINE COORDINATES OF VERTICES AT TOP
+    }
+    std::vector<int> vertexIDs(numberOfVertices);
+    
+    participantSolid.setMeshVertices(meshName, coords, vertexIDs);
+
+    // ----- initialize read- and write-data
+    std::vector<double> temperature(numberOfVertices * meshDim);    // read-data
+    std::vector<double> 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<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_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<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;
+
+
+    // ----- free data, close communication -----
+    participantSolid.finalize();
+
+
+
+    return 0;
+}

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