max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-15 18:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'hannes-philipp' of git.gfz-potsdam.de:naaice/tug into hannes-philipp

Browse Source
This commit is contained in:
philippun 2023-08-22 13:32:13 +02:00
commit 27829a1463
29 changed files with 8742 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
.readthedocs.yaml Normal file
View File

@ -0,0 +1,32 @@
# .readthedocs.yaml
# Read the Docs configuration file
# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
# Required
version: 2
# Set the OS, Python version and other tools you might need
build:
os: ubuntu-22.04
tools:
python: "3.11"
# You can also specify other tool versions:
# nodejs: "19"
# rust: "1.64"
# golang: "1.19"
# Build documentation in the "docs/" directory with Sphinx
sphinx:
configuration: docs_sphinx/conf.py
# Optionally build your docs in additional formats such as PDF and ePub
formats:
- pdf
# - epub
# Optional but recommended, declare the Python requirements required
# to build your documentation
# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
python:
install:
- requirements: docs_sphinx/requirements.txt

4
docs_sphinx/Boundary.rst Normal file
View File

@ -0,0 +1,4 @@
Boundary
========
.. doxygenclass:: Boundary

2771
docs_sphinx/Doxyfile.in Normal file
View File

File diff suppressed because it is too large Load Diff

4
docs_sphinx/Grid.rst Normal file
View File

@ -0,0 +1,4 @@
Grid
====
.. doxygenclass:: Grid

20
docs_sphinx/Makefile Normal file
View File

@ -0,0 +1,20 @@
# Minimal makefile for Sphinx documentation
#
# You can set these variables from the command line, and also
# from the environment for the first two.
SPHINXOPTS ?=
SPHINXBUILD ?= sphinx-build
SOURCEDIR = .
BUILDDIR = _build
# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
.PHONY: help Makefile
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

4
docs_sphinx/Simulation.rst Normal file
View File

@ -0,0 +1,4 @@
Simulation
==========
.. doxygenclass:: Simulation

4
docs_sphinx/_build/html/_images/activity_diagram_run.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 15 KiB

4
docs_sphinx/_build/html/_images/class_diagram.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 23 KiB

224
docs_sphinx/_build/vis.html Normal file
View File

@ -0,0 +1,224 @@
<!DOCTYPE html>
<div id="container"></div>
<form>
<label for="c_file">Select an output file: </label>
<input type="file" id="c_file" name="c_file" accept="text/csv" />
</form>
<script type="module">
import * as d3 from "https://cdn.jsdelivr.net/npm/d3@7/+esm";
// Declare the chart dimensions and margins.
const width = 600;
const height = 600;
const svgMargin = { top: 50, right: 50, bottom: 50, left: 50 };
var gridWidth = width - svgMargin.left - svgMargin.right;
var gridHeight = height - svgMargin.top - svgMargin.bottom;
const gap = 1;
var state = 0;
var numIterations = 0;
var cellXcount = 0;
var cellYcount = 0;
var concentrations = [];
// Create the outer SVG container.
const svg = d3.create("svg")
.attr("width", width)
.attr("height", height);
svg.append("rect")
.attr("x", 0)
.attr("y", 0)
.attr("width", width)
.attr("height", height)
.attr("fill-opacity", 0)
.attr("stroke-width", 4)
.attr("stroke", "black");
var grid = svg.append("g")
.attr("class", "grid")
.attr("transform", `translate(${svgMargin.left}, ${svgMargin.top})`);
// color scale
var colorScale = d3.scaleLinear()
.range(["#d0d0ec", "#7f0000"])
.domain([1, 2000]);
function getColor(c, max) {
if (c == -1) {
return "#000000";
}
colorScale.domain([0, max]);
return colorScale(c);
}
function calcMaxConcentrationInIteration(state) {
var maxRow = concentrations[state].map(function(row){ return Math.max.apply(Math, row); });
var maxC = Math.max.apply(null, maxRow);
return maxC;
}
// Load data from file
async function createVisualFromFile(csv_input) {
// console.log(csv_input);
// var data = await d3.text(csv_input).then(function(csv){
// return d3.dsvFormat(" ").parseRows(csv);
// })
var data = d3.dsvFormat(" ").parseRows(csv_input);
// console.log(data);
var leftBoundary = data[0];
var rightBoundary = data[1];
var topBoundary = data[2];
var bottomBoundary = data[3];
cellXcount = topBoundary.length;
cellYcount = leftBoundary.length;
concentrations = []; // reset concentrations
grid.selectAll("rect").remove();
svg.selectAll("rect").remove();
state = 0;
console.log(topBoundary);
var cellWidth = gridWidth / cellXcount;
var cellHeight = gridHeight / cellYcount;
var iteration = [];
numIterations = (data.length - 6) / (cellYcount + 2);
console.log(numIterations);
for (let i = 0; i < numIterations; i++) {
iteration = [];
for (let j = 0; j < cellYcount; j++) {
iteration.push(data[i * (cellYcount + 2) + 6 + j])
}
concentrations.push(iteration);
}
console.log(concentrations);
var maxC = calcMaxConcentrationInIteration(state);
// Create grid
for (let i = 0; i < cellYcount; i++) {
for (let j = 0; j < cellXcount; j++) {
grid.append("rect")
.attr("x", i * cellWidth + gap/2)
.attr("y", j * cellHeight + gap/2)
.attr("width", cellWidth - gap)
.attr("height", cellHeight - gap)
.attr("fill", getColor(concentrations[state][i][j], maxC));
}
}
// Create Boundaries
// left and right
for (let j = 0; j < cellYcount; j++) {
svg.append("rect")
.attr("x", svgMargin.left - 10)
.attr("y", svgMargin.top + j * cellHeight + gap/2)
.attr("width", 7)
.attr("height", cellHeight - gap)
.attr("fill", getColor(leftBoundary[j], maxC));
svg.append("rect")
.attr("x", width - svgMargin.right + 3)
.attr("y", svgMargin.top + j * cellHeight + gap/2)
.attr("width", 7)
.attr("height", cellHeight - gap)
.attr("fill", getColor(rightBoundary[j], maxC));
}
// top and bottom
for (let i = 0; i < cellXcount; i++) {
svg.append("rect")
.attr("x", svgMargin.left + i * cellWidth + gap/2)
.attr("y", svgMargin.top - 10)
.attr("width", cellWidth - gap)
.attr("height", 7)
.attr("fill", getColor(topBoundary[i], maxC));
svg.append("rect")
.attr("x", svgMargin.left + i * cellWidth + gap/2)
.attr("y", height - svgMargin.bottom + 3)
.attr("width", cellWidth - gap)
.attr("height", 7)
.attr("fill", getColor(bottomBoundary[i], maxC));
}
}
function updateGrid(new_state) {
var maxC = calcMaxConcentrationInIteration(new_state);
console.log(maxC);
grid.selectAll("rect")
.attr("fill", function (d,i) {
var row = Math.floor(i/20);
var col = i%20;
return getColor(concentrations[new_state][row][col], maxC);
})
}
// key events for changing grid iteration state
addEventListener("keydown", (event) => {
if (event.isComposing || event.keyCode === 81) {
if (state > 0) {
state -= 1;
updateGrid(state);
}
}
if (event.isComposing || event.keyCode === 69) {
if (state < numIterations-1) {
state += 1;
updateGrid(state);
}
}
if (event.isComposing || event.keyCode === 65) {
if (state > 9) {
state -= 10;
} else {
state = 0;
}
updateGrid(state);
}
if (event.isComposing || event.keyCode === 68) {
if (state < numIterations-10) {
state += 10;
} else {
state = numIterations-1;
}
updateGrid(state);
}
});
async function logFile() {
var file = document.getElementById("c_file").files[0];
const reader = new FileReader();
reader.readAsText(file);
console.log(reader.result);
}
function reader(file, callback) {
const fr = new FileReader();
fr.onload = () => callback(null, fr.result);
fr.onerror = (err) => callback(err);
fr.readAsText(file);
}
// document.getElementById("c_file").addEventListener("change", logFile);
document.getElementById("c_file").addEventListener("change", (evt) => {
// No files, do nothing.
if (!evt.target.files) {
return;
}
reader(evt.target.files[0], (err, res) => {
console.log(res); // Base64 `data:image/...` String result.
createVisualFromFile(res);
});
});
// Append the SVG element.
container.append(svg.node());
</script>

97
docs_sphinx/conf.py Normal file
View File

@ -0,0 +1,97 @@
# Configuration file for the Sphinx documentation builder.
#
# This file only contains a selection of the most common options. For a full
# list see the documentation:
# https://www.sphinx-doc.org/en/master/usage/configuration.html
# -- Path setup --------------------------------------------------------------
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
# import os
# import sys
# sys.path.insert(0, os.path.abspath('.'))
from sphinx.builders.html import StandaloneHTMLBuilder
import subprocess, os
# Doxygen
subprocess.call('doxygen Doxyfile.in', shell=True)
# -- Project information -----------------------------------------------------
project = 'TUG'
copyright = 'MIT'
author = 'Philipp Ungrund, Hannes Signer'
# -- General configuration ---------------------------------------------------
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = [
'sphinx.ext.autodoc',
'sphinx.ext.intersphinx',
'sphinx.ext.autosectionlabel',
'sphinx.ext.todo',
'sphinx.ext.coverage',
'sphinx.ext.mathjax',
'sphinx.ext.ifconfig',
'sphinx.ext.viewcode',
'sphinx.ext.inheritance_diagram',
'breathe'
]
html_baseurl = "/index.html"
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This pattern also affects html_static_path and html_extra_path.
exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
highlight_language = 'c++'
# -- Options for HTML output -------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
html_theme = 'sphinx_rtd_theme'
html_theme_options = {
'canonical_url': '',
'analytics_id': '', # Provided by Google in your dashboard
'display_version': True,
'prev_next_buttons_location': 'bottom',
'style_external_links': False,
'logo_only': False,
# Toc options
'collapse_navigation': True,
'sticky_navigation': True,
'navigation_depth': 4,
'includehidden': True,
'titles_only': False
}
# html_logo = ''
# github_url = ''
# html_baseurl = ''
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
html_logo = "images/tug_logo.svg"
# -- Breathe configuration -------------------------------------------------
breathe_projects = {
"Tug": "_build/xml/"
}
breathe_default_project = "Tug"
breathe_default_members = ('members', 'undoc-members')

38
docs_sphinx/developer.rst Normal file
View File

@ -0,0 +1,38 @@
Developer Guide
===============
=========================
Class Diagram of user API
=========================
The following graphic shows the class diagram of the user API. The FTCS and
BTCS functionalities are externally outsourced and not visible to the user.
.. image:: images/class_diagram.svg
:width: 2000
:alt: Class diagram for the user API
====================================================
Activity Diagram for run routine in simulation class
====================================================
The following activity diagram represents the actions when the run method is called within the simulation class.
For better distinction, the activities of the calculation methods FTCS and BTCS are shown in two separate activity diagrams.
.. image:: images/activity_diagram_run.svg
:width: 2000
:alt: Activity diagram for the run method in the simulation class
**Activity Diagram for FTCS method**
.. image:: images/activity_diagram_FTCS.svg
:width: 400
:alt: Activity diagram for the FTCS method
**Activity Diagram for BTCS method**
.. image:: images/activity_diagram_BTCS.svg
:width: 400
:alt: Activity diagram for the BTCS method

2
docs_sphinx/developper.rst Normal file
View File

@ -0,0 +1,2 @@
Developper API
==============

62
docs_sphinx/examples.rst Normal file
View File

@ -0,0 +1,62 @@
Examples
========
At this point, some typical commented examples are presented to illustrate how Tug works.
In general, each simulation is divided into three blocks:
- the initialization of the grid, which is to be simulated
- the setting of the boundary conditions
- the setting of the simulation parameters and the start of the simulation
Two dimensional grid with constant boundaries and FTCS method
-------------------------------------------------------------
**Initialization of the grid**
For example, the initalization of a grid with 20 by 20 cells and a domain size (physical extent of the grid) of
also 20 by 20 length units can be done as follows. The setting of the domain is optional here and is set to the
same size as the number of cells in the standard case. As seen in the code, the cells of the grid are set to an
initial value of 0 and only in the upper left corner (0,0) the starting concentration is set to the value 20.
.. code-block:: cpp
int row = 20
int col = 20;
Grid grid = Grid(row,col);
grid.setDomain(row, col);
MatrixXd concentrations = MatrixXd::Constant(row,col,0);
concentrations(0,0) = 20;
grid.setConcentrations(concentrations);
**Setting of the boundary conditions**
First, a boundary class is created and then the corresponding boundary conditions are set. In this case, all four sides
of the grid are set as constant edges with a concentration of 0.
.. code-block:: cpp
Boundary bc = Boundary(grid);
bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
**Setting of the simulation parameters and simulation start**
In the last block, a simulation class is created and the objects of the grid and the boundary conditions are passed. The solution
method is also specified (either FCTS or BTCS). Furthermore, the desired time step and the number of iterations are set. The penultimate
parameter specifies the output of the simulated results in a CSV file. In the present case, the result of each iteration step is written
one below the other into the corresponding CSV file.
.. code-block:: cpp
Simulation simulation = Simulation(grid, bc, FTCS_APPROACH);
simulation.setTimestep(0.1);
simulation.setIterations(1000);
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
simulation.run();
Setting special boundary conditions on individual cells
-------------------------------------------------------

4
docs_sphinx/images/activity_diagram_BTCS.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 35 KiB

4
docs_sphinx/images/activity_diagram_FTCS.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 19 KiB

4
docs_sphinx/images/activity_diagram_run.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 15 KiB

4
docs_sphinx/images/class_diagram.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 23 KiB

2629
docs_sphinx/images/tug_logo.svg Normal file
View File

File diff suppressed because it is too large Load Diff
Side by Side

After

Width:  |  Height:  |  Size: 138 KiB

37
docs_sphinx/index.rst Normal file
View File

@ -0,0 +1,37 @@
.. Tug documentation master file, created by
sphinx-quickstart on Mon Aug 14 11:30:23 2023.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
Welcome to Tug's documentation!
===============================
Welcome to the documentation of the TUG project, a simulation program
for solving one- and two-dimensional diffusion problems with heterogeneous diffusion coefficients, more
generally, for solving the following differential equation
.. math::
\frac{\partial C}{\partial t} = \alpha_x \frac{\partial^2 C}{\partial x^2} + \alpha_y \frac{\partial^2 C}{\partial y^2}.
.. toctree::
:maxdepth: 2
:caption: Contents:
Indices and tables
==================
* :ref:`genindex`
* :ref:`search`
Table of Contents
^^^^^^^^^^^^^^^^^
.. toctree::
:maxdepth: 2
self
installation
theory
user
developer
examples
visualization

3
docs_sphinx/installation.rst Normal file
View File

@ -0,0 +1,3 @@
Installation
============

2
docs_sphinx/requirements.txt Normal file
View File

@ -0,0 +1,2 @@
breathe
sphinx-rtd-theme

15
docs_sphinx/theory.rst Normal file
View File

@ -0,0 +1,15 @@
Theoretical Foundations
=======================
=====================
The Diffusion Problem
=====================
================
Numerical Solver
================
**Backward Time-Centered Space (BTCS) Method**
**Forward Time-Centered Space (BTCS) Method**

2629
docs_sphinx/tug_logo.svg Normal file
View File

File diff suppressed because it is too large Load Diff
Side by Side

After

Width:  |  Height:  |  Size: 138 KiB

9
docs_sphinx/user.rst Normal file
View File

@ -0,0 +1,9 @@
User API
========
.. toctree::
:maxdepth: 2
Grid
Boundary
Simulation

3
docs_sphinx/visualization.rst Normal file
View File

@ -0,0 +1,3 @@
Visualization
=============

50
examples/profiling.cpp Normal file
View File

@ -0,0 +1,50 @@
#include <tug/Simulation.hpp>
#include <iostream>
#include <fstream>
int main(int argc, char *argv[]) {
int n[4] = {10, 20, 50, 100};
int iterations[1] = {100};
int repetition = 10;
ofstream myfile;
myfile.open("btcs_time_measurement_openmp_10.csv");
for (int i = 0; i < size(n); i++){
cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
myfile << "Grid size: " << n[i] << " x " << n[i] << endl;
for(int j = 0; j < size(iterations); j++){
cout << "Iterations: " << iterations[j] << endl;
for (int k = 0; k < repetition; k++){
cout << "Wiederholung: " << k << endl;
Grid grid = Grid(n[i], n[i]);
grid.setDomain(n[i], n[i]);
MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);
concentrations(5,5) = 1;
grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
Boundary bc = Boundary(grid);
Simulation sim = Simulation(grid, bc, BTCS_APPROACH);
sim.setTimestep(0.001);
sim.setIterations(iterations[j]);
sim.setOutputCSV(CSV_OUTPUT_OFF);
auto begin = std::chrono::high_resolution_clock::now();
sim.run();
auto end = std::chrono::high_resolution_clock::now();
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
myfile << milliseconds.count() << endl;
}
}
cout << endl;
myfile << endl;
}
myfile.close();
}

51
examples/profiling_openmp.cpp Normal file
View File

@ -0,0 +1,51 @@
#include <tug/Simulation.hpp>
#include <iostream>
#include <fstream>
int main(int argc, char *argv[]) {
int n[4] = {10, 20, 50, 100};
int iterations[1] = {100};
int repetition = 10;
ofstream myfile;
myfile.open("time_measure_experiment_openmp_thread_6.csv");
for (int i = 0; i < size(n); i++){
cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
for(int j = 0; j < size(iterations); j++){
cout << "Iterations: " << iterations[j] << endl;
myfile << "Iterations: " << iterations[j] << endl;
for (int k = 0; k < repetition; k++){
cout << "Wiederholung: " << k << endl;
Grid grid = Grid(n[i], n[i]);
grid.setDomain(n[i], n[i]);
MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);
concentrations(5,5) = 1;
grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
Boundary bc = Boundary(grid);
Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
sim.setTimestep(0.001);
sim.setIterations(iterations[j]);
sim.setOutputCSV(CSV_OUTPUT_OFF);
auto begin = std::chrono::high_resolution_clock::now();
sim.run();
auto end = std::chrono::high_resolution_clock::now();
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
myfile << milliseconds.count() << endl;
}
}
cout << endl;
myfile << endl;
}
myfile.close();
}

36
src/BTCSv2.cpp
View File

@ -7,6 +7,9 @@
#include "FTCS.cpp" #include "FTCS.cpp"
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#include <omp.h>
#define NUM_THREADS_BTCS 1
using namespace Eigen; using namespace Eigen;
@ -62,7 +65,7 @@ static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &
// creates coefficient matrix for next time step from alphas in x-direction // creates coefficient matrix for next time step from alphas in x-direction
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> bcLeft, vector<BoundaryElement> bcRight, int numCols, int rowIndex, double sx) { static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, int &numCols, int &rowIndex, double &sx) {
// square matrix of column^2 dimension for the coefficients // square matrix of column^2 dimension for the coefficients
SparseMatrix<double> cm(numCols, numCols); SparseMatrix<double> cm(numCols, numCols);
@ -84,6 +87,7 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryEl
// inner columns // inner columns
int n = numCols-1; int n = numCols-1;
#pragma omp parallel for num_threads(NUM_THREADS_BTCS)
for (int i = 1; i < n; i++) { for (int i = 1; i < n; i++) {
cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)); cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i));
cm.insert(i,i) = 1 + sx * ( cm.insert(i,i) = 1 + sx * (
@ -190,7 +194,7 @@ static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentra
static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight,
vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom, vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom,
int length, int rowIndex, double sx, double sy) { int &length, int &rowIndex, double &sx, double &sy) {
VectorXd sv(length); VectorXd sv(length);
int numRows = concentrations.rows(); int numRows = concentrations.rows();
@ -198,6 +202,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
// inner rows // inner rows
if (rowIndex > 0 && rowIndex < numRows-1) { if (rowIndex > 0 && rowIndex < numRows-1) {
#pragma omp parallel for num_threads(NUM_THREADS_BTCS)
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
* concentrations(rowIndex+1,i) * concentrations(rowIndex+1,i)
@ -215,6 +220,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
// first row // first row
if (rowIndex == 0) { if (rowIndex == 0) {
#pragma omp parallel for num_threads(NUM_THREADS_BTCS)
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
type = bcTop[i].getType(); type = bcTop[i].getType();
if (type == BC_TYPE_CONSTANT) { if (type == BC_TYPE_CONSTANT) {
@ -229,6 +235,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
// last row // last row
if (rowIndex == numRows-1) { if (rowIndex == numRows-1) {
#pragma omp parallel for num_threads(NUM_THREADS_BTCS)
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
type = bcBottom[i].getType(); type = bcBottom[i].getType();
if (type == BC_TYPE_CONSTANT) { if (type == BC_TYPE_CONSTANT) {
@ -258,6 +265,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
// solver for linear equation system; A corresponds to coefficient matrix, // solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector // b to the solution vector
static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) { static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
SparseLU<SparseMatrix<double>> solver; SparseLU<SparseMatrix<double>> solver;
solver.analyzePattern(A); solver.analyzePattern(A);
solver.factorize(A); solver.factorize(A);
@ -281,7 +289,8 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
MatrixXd concentrations = grid.getConcentrations(); MatrixXd concentrations = grid.getConcentrations();
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, 0, sx); // this is exactly same as in 2D int rowIndex = 0;
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
b(i) = concentrations(0,i); b(i) = concentrations(0,i);
} }
@ -323,34 +332,39 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
MatrixXd concentrations = grid.getConcentrations(); MatrixXd concentrations = grid.getConcentrations();
#pragma omp parallel for num_threads(NUM_THREADS_BTCS) private(A, b, row_t1)
for (int i = 0; i < rowMax; i++) { for (int i = 0; i < rowMax; i++) {
A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
bcTop, bcBottom, colMax, i, sx, sy); bcTop, bcBottom, colMax, i, sx, sy);
SparseLU<SparseMatrix<double>> solver;
row_t1 = solve(A, b); row_t1 = solve(A, b);
for (int j = 0; j < colMax; j++) { concentrations_t1.row(i) = row_t1;
concentrations_t1(i,j) = row_t1(j); // can potentially be improved by using Eigen method
}
} }
concentrations_t1.transposeInPlace(); concentrations_t1.transposeInPlace();
concentrations.transposeInPlace(); concentrations.transposeInPlace();
alphaX.transposeInPlace(); alphaX.transposeInPlace();
alphaY.transposeInPlace(); alphaY.transposeInPlace();
#pragma omp parallel for num_threads(NUM_THREADS_BTCS) private(A, b, row_t1)
for (int i = 0; i < colMax; i++) { for (int i = 0; i < colMax; i++) {
// swap alphas, boundary conditions and sx/sy for column-wise calculation // swap alphas, boundary conditions and sx/sy for column-wise calculation
A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy);
b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
bcLeft, bcRight, rowMax, i, sy, sx); bcLeft, bcRight, rowMax, i, sy, sx);
row_t1 = solve(A, b); row_t1 = solve(A, b);
for (int j = 0; j < rowMax; j++) { concentrations.row(i) = row_t1;
concentrations(i,j) = row_t1(j); // can potentially be improved by using Eigen method
}
} }
concentrations.transposeInPlace(); concentrations.transposeInPlace();
grid.setConcentrations(concentrations); grid.setConcentrations(concentrations);

12
src/FTCS.cpp
View File

@ -11,6 +11,8 @@
#include <iostream> #include <iostream>
#include <omp.h> #include <omp.h>
#define NUM_THREADS 6
using namespace std; using namespace std;
@ -269,7 +271,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// inner cells // inner cells
// these are independent of the boundary condition type // these are independent of the boundary condition type
// omp_set_num_threads(10); // omp_set_num_threads(10);
#pragma omp parallel for #pragma omp parallel for num_threads(NUM_THREADS)
for (int row = 1; row < rowMax-1; row++) { for (int row = 1; row < rowMax-1; row++) {
for (int col = 1; col < colMax-1; col++) { for (int col = 1; col < colMax-1; col++) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col) concentrations_t1(row, col) = grid.getConcentrations()(row, col)
@ -289,7 +291,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// left without corners / looping over rows // left without corners / looping over rows
// hold column constant at index 0 // hold column constant at index 0
int col = 0; int col = 0;
#pragma omp parallel for #pragma omp parallel for num_threads(NUM_THREADS)
for (int row = 1; row < rowMax-1; row++) { for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row, col) = grid.getConcentrations()(row,col) concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol) + timestep / (deltaCol*deltaCol)
@ -306,7 +308,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// right without corners / looping over rows // right without corners / looping over rows
// hold column constant at max index // hold column constant at max index
col = colMax-1; col = colMax-1;
#pragma omp parallel for #pragma omp parallel for num_threads(NUM_THREADS)
for (int row = 1; row < rowMax-1; row++) { for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col) concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol) + timestep / (deltaCol*deltaCol)
@ -324,7 +326,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// top without corners / looping over columns // top without corners / looping over columns
// hold row constant at index 0 // hold row constant at index 0
int row = 0; int row = 0;
#pragma omp parallel for #pragma omp parallel for num_threads(NUM_THREADS)
for (int col=1; col<colMax-1;col++){ for (int col=1; col<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col) concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow) + timestep / (deltaRow*deltaRow)
@ -341,7 +343,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// bottom without corners / looping over columns // bottom without corners / looping over columns
// hold row constant at max index // hold row constant at max index
row = rowMax-1; row = rowMax-1;
#pragma omp parallel for #pragma omp parallel for num_threads(NUM_THREADS)
for(int col=1; col<colMax-1;col++){ for(int col=1; col<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col) concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow) + timestep / (deltaRow*deltaRow)