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

This commit is contained in:
Marco De Lucia 2023-08-24 17:29:58 +02:00
commit 2294922a3e
3 changed files with 561 additions and 18 deletions

View File

@ -8,8 +8,8 @@ int main(int argc, char *argv[]) {
// **************
// profiler::startListen();
// create a grid with a 20 x 20 field
int row = 20;
int col = 20;
int row = 40;
int col = 50;
Grid grid = Grid(row,col);
// (optional) set the domain, e.g.:
@ -19,7 +19,7 @@ int main(int argc, char *argv[]) {
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value
// grid.setConcentrations(concentrations);
MatrixXd concentrations = MatrixXd::Constant(row,col,0);
concentrations(0,0) = 2000;
concentrations(10,10) = 2000;
grid.setConcentrations(concentrations);
@ -67,7 +67,7 @@ int main(int argc, char *argv[]) {
simulation.setTimestep(0.1); // timestep
// set the number of iterations
simulation.setIterations(100);
simulation.setIterations(300);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_XTREME);

View File

@ -15,7 +15,7 @@ using namespace Eigen;
// calculates coefficient for left boundary in constant case
static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, double &sx) {
static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff;
double rightCoeff;
@ -28,7 +28,7 @@ static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int
// calculates coefficient for left boundary in closed case
static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, double &sx) {
static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff;
double rightCoeff;
@ -40,7 +40,7 @@ static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &r
// calculates coefficient for right boundary in constant case
static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, int n, double sx) {
double leftCoeff;
double centerCoeff;
@ -53,7 +53,7 @@ static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int
// calculates coefficient for right boundary in closed case
static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
double leftCoeff;
double centerCoeff;
@ -65,7 +65,7 @@ static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &
// 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
SparseMatrix<double> cm(numCols, numCols);
@ -119,7 +119,7 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryEl
// calculates explicity concentration at top boundary in constant case
static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations,
MatrixXd &alpha, vector<BoundaryElement> &bcTop, int &rowIndex, int &i, double &sy) {
MatrixXd &alpha, vector<BoundaryElement> &bcTop, int rowIndex, int i, double sy) {
double c;
c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
@ -138,7 +138,7 @@ static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrat
// calculates explicit concentration at top boundary in closed case
static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations,
MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
MatrixXd &alpha, int rowIndex, int i, double sy) {
double c;
c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
@ -155,7 +155,7 @@ static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentratio
// calculates explicit concentration at bottom boundary in constant case
static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations,
MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int &rowIndex, int &i, double &sy) {
MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
double c;
c = sy * alpha(rowIndex,i) * bcBottom[i].getValue()
@ -174,7 +174,7 @@ static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concent
// calculates explicit concentration at bottom boundary in closed case
static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations,
MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
MatrixXd &alpha, int rowIndex, int i, double sy) {
double c;
c = (
@ -193,7 +193,7 @@ static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentra
static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight,
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);
int numRows = concentrations.rows();
@ -316,7 +316,7 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
// BTCS solution for 1D grid
static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b)) {
static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b)) {
int length = grid.getLength();
double sx = timestep / (grid.getDelta() * grid.getDelta());
@ -353,7 +353,7 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solve
// BTCS solution for 2D grid
static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b), int &numThreads) {
static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b), int numThreads) {
int rowMax = grid.getRow();
int colMax = grid.getCol();
double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
@ -413,7 +413,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solve
// entry point for EigenLU solver; differentiate between 1D and 2D grid
static void BTCS_LU(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
} else if (grid.getDim() == 2) {
@ -424,7 +424,7 @@ static void BTCS_LU(Grid &grid, Boundary &bc, double &timestep, int &numThreads)
}
// entry point for Thomas algorithm solver; differentiate 1D and 2D grid
static void BTCS_Thomas(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
} else if (grid.getDim() == 2) {

View File

@ -0,0 +1,543 @@
<!DOCTYPE html>
<html lang="en-US">
<head>
<!-- <link rel="stylesheet" href="style.css"> -->
<style>
body {
background-color: white;
}
:root {
--max-width: 1440px;
}
.main-wrapper {
display: grid;
grid-template-areas: "sidebar main";
grid-template-columns: minmax(0,1fr) minmax(0,2.5fr);
grid-column-gap: 50px;
margin: 0 auto;
max-width: var(--max-width);
}
#state-info {
background-color: lightgray;
width: 120px;
text-align: center;
border: 2px solid black;
border-radius: 5px;
font-family: Arial, Helvetica, sans-serif;
font-size: 25px;
margin-bottom: 20px;
}
</style>
</head>
<body>
<div id="root">
<div class="main-wrapper">
<div class="sidebar-container">
<aside id="sidebar">
<div class="menu">
<form>
<label for="c_file">File (XTREME-format): </label>
<input type="file" id="c_file" name="c_file" accept="text/csv" />
</form>
</div>
<div class="info">
<p>
Use the keys <i>q</i> and <i>e</i> for increasing and decreasing the state by <b>1</b>.
The keys <i>a</i> and <i>d</i> for increasing and decresing the state by <b>10</b>.
</p>
<div id="state-info">
State: -
</div>
<div id="legend-info">
</div>
</div>
</aside>
</div>
<div class="content-container">
<main id="content" class="main-content">
</main>
</div>
</div>
</div>
</body>
</html>
<script type="module">
import * as d3 from "https://cdn.jsdelivr.net/npm/d3@7/+esm";
// layout information; constants
const width = document.getElementById('content').offsetWidth;
const height = document.getElementById('content').offsetWidth;
const svgMargin = { top: 10, right: 10, bottom: 10, left: 10 };
const boundary_gap = 1;
const gridWidth = width - svgMargin.left - svgMargin.right - 2*boundary_gap;
const gridHeight = height - svgMargin.top - svgMargin.bottom - 2*boundary_gap;
// simulation information
var state = 0; // iteration state
var iteration_count = 0;
var number_rows = 0;
var number_cols = 0;
var max_concentration = 0;
var data = [];
var leftb, rightb, topb, bottomb;
// test data; to be deleted later on
var test_data = [{"row": 0, "col": 0, "value": 1},
{"row": 0, "col": 1, "value": 1},
{"row": 0, "col": 2, "value": 1},
{"row": 0, "col": 3, "value": 1},
{"row": 1, "col": 0, "value": 1},
{"row": 1, "col": 1, "value": 1},
{"row": 1, "col": 2, "value": 1},
{"row": 1, "col": 3, "value": 1},
{"row": 2, "col": 0, "value": 1},
{"row": 2, "col": 1, "value": 1},
{"row": 2, "col": 2, "value": 1},
{"row": 2, "col": 3, "value": 1},
{"row": 3, "col": 0, "value": 1},
{"row": 3, "col": 1, "value": 1},
{"row": 3, "col": 2, "value": 1},
{"row": 3, "col": 3, "value": 1}];
//////////////////////////////////////////
////////
//////////////////////////////////////////
// Copyright 2021, Observable Inc.
// Released under the ISC license.
// https://observablehq.com/@d3/color-legend
function Legend(color, {
title,
tickSize = 6,
width = 320,
height = 44 + tickSize,
marginTop = 18,
marginRight = 0,
marginBottom = 16 + tickSize,
marginLeft = 0,
ticks = width / 64,
tickFormat,
tickValues
} = {}) {
function ramp(color, n = 256) {
const canvas = document.createElement("canvas");
canvas.width = n;
canvas.height = 1;
const context = canvas.getContext("2d");
for (let i = 0; i < n; ++i) {
context.fillStyle = color(i / (n - 1));
context.fillRect(i, 0, 1, 1);
}
return canvas;
}
const svg = d3.create("svg")
.attr("width", width)
.attr("height", height)
.attr("viewBox", [0, 0, width, height])
.style("overflow", "visible")
.style("display", "block");
let tickAdjust = g => g.selectAll(".tick line").attr("y1", marginTop + marginBottom - height);
let x;
// Continuous
if (color.interpolate) {
const n = Math.min(color.domain().length, color.range().length);
x = color.copy().rangeRound(d3.quantize(d3.interpolate(marginLeft, width - marginRight), n));
svg.append("image")
.attr("x", marginLeft)
.attr("y", marginTop)
.attr("width", width - marginLeft - marginRight)
.attr("height", height - marginTop - marginBottom)
.attr("preserveAspectRatio", "none")
.attr("xlink:href", ramp(color.copy().domain(d3.quantize(d3.interpolate(0, 1), n))).toDataURL());
}
// Sequential
else if (color.interpolator) {
x = Object.assign(color.copy()
.interpolator(d3.interpolateRound(marginLeft, width - marginRight)),
{range() { return [marginLeft, width - marginRight]; }});
svg.append("image")
.attr("x", marginLeft)
.attr("y", marginTop)
.attr("width", width - marginLeft - marginRight)
.attr("height", height - marginTop - marginBottom)
.attr("preserveAspectRatio", "none")
.attr("xlink:href", ramp(color.interpolator()).toDataURL());
// scaleSequentialQuantile doesnt implement ticks or tickFormat.
if (!x.ticks) {
if (tickValues === undefined) {
const n = Math.round(ticks + 1);
tickValues = d3.range(n).map(i => d3.quantile(color.domain(), i / (n - 1)));
}
if (typeof tickFormat !== "function") {
tickFormat = d3.format(tickFormat === undefined ? ",f" : tickFormat);
}
}
}
// Threshold
else if (color.invertExtent) {
const thresholds
= color.thresholds ? color.thresholds() // scaleQuantize
: color.quantiles ? color.quantiles() // scaleQuantile
: color.domain(); // scaleThreshold
const thresholdFormat
= tickFormat === undefined ? d => d
: typeof tickFormat === "string" ? d3.format(tickFormat)
: tickFormat;
x = d3.scaleLinear()
.domain([-1, color.range().length - 1])
.rangeRound([marginLeft, width - marginRight]);
svg.append("g")
.selectAll("rect")
.data(color.range())
.join("rect")
.attr("x", (d, i) => x(i - 1))
.attr("y", marginTop)
.attr("width", (d, i) => x(i) - x(i - 1))
.attr("height", height - marginTop - marginBottom)
.attr("fill", d => d);
tickValues = d3.range(thresholds.length);
tickFormat = i => thresholdFormat(thresholds[i], i);
}
// Ordinal
else {
x = d3.scaleBand()
.domain(color.domain())
.rangeRound([marginLeft, width - marginRight]);
svg.append("g")
.selectAll("rect")
.data(color.domain())
.join("rect")
.attr("x", x)
.attr("y", marginTop)
.attr("width", Math.max(0, x.bandwidth() - 1))
.attr("height", height - marginTop - marginBottom)
.attr("fill", color);
tickAdjust = () => {};
}
svg.append("g")
.attr("transform", `translate(0,${height - marginBottom})`)
.call(d3.axisBottom(x)
.ticks(ticks, typeof tickFormat === "string" ? tickFormat : undefined)
.tickFormat(typeof tickFormat === "function" ? tickFormat : undefined)
.tickSize(tickSize)
.tickValues(tickValues))
.call(tickAdjust)
.call(g => g.select(".domain").remove())
.call(g => g.append("text")
.attr("x", marginLeft)
.attr("y", marginTop + marginBottom - height - 6)
.attr("fill", "currentColor")
.attr("text-anchor", "start")
.attr("font-weight", "bold")
.attr("class", "title")
.text(title));
return svg.node();
}
//////////////////////////////////////////
////////
//////////////////////////////////////////
// max concentration
function calc_max_concentration(state) {
// var maxRow = data[state].map(function(row){ return Math.max.apply(Math, row); });
// var maxC = Math.max.apply(null, maxRow);
var maxC = Math.max(...data[state].map(d => d.value))
return maxC;
}
// scales
var col, row, color;
function create_n_update_scales() {
col = d3.scaleBand()
.range([0, gridWidth])
.domain([...Array(number_cols).keys()])
.padding(0.05);
row = d3.scaleBand()
.range([0, gridHeight])
.domain([...Array(number_rows).keys()])
.padding(0.05);
color = d3.scaleLinear()
.range(["#d0d0ec", "#7f0000"])
.domain([0, max_concentration]);
}
// grid
function draw_n_update_grid(grid_data) {
document.getElementById('state-info').innerHTML = `State: ${state}`;
// Grid rects
grid.selectAll()
.data(grid_data, (d,i) => d.row+':'+d.col)
.enter()
.append("rect")
.attr("x", (d,i) => col(d.col))
.attr("y", (d,i) => row(d.row))
.attr("width", col.bandwidth())
.attr("height", row.bandwidth())
.attr("fill", (d,i) => color(d.value))
.on("mouseover", mouseover)
.on("mousemove", mousemove)
.on("mouseleave", mouseleave);
}
// boundaries
function draw_boundaries() {
left_side.selectAll("rect")
.data(leftb)
.join("rect")
.attr("x", (d,i) => 0)
.attr("y", (d,i) => row(i))
.attr("width", svgMargin.left)
.attr("height", row.bandwidth())
.attr("fill", function(d,i) {
if (d == -1) {
return "#000000";
} else {
return color(d);
}
});
right_side.selectAll("rect")
.data(rightb)
.join("rect")
.attr("x", (d,i) => 0)
.attr("y", (d,i) => row(i))
.attr("width", svgMargin.right)
.attr("height", row.bandwidth())
.attr("fill", function(d,i) {
if (d == -1) {
return "#000000";
} else {
return color(d);
}
});
top_side.selectAll("rect")
.data(topb)
.join("rect")
.attr("x", (d,i) => col(i))
.attr("y", (d,i) => 0)
.attr("width", col.bandwidth())
.attr("height", svgMargin.top)
.attr("fill", function(d,i) {
if (d == -1) {
return "#000000";
} else {
return color(d);
}
});
bottom_side.selectAll("rect")
.data(bottomb)
.join("rect")
.attr("x", (d,i) => col(i))
.attr("y", (d,i) => 0)
.attr("width", col.bandwidth())
.attr("height", svgMargin.bottom)
.attr("fill", function(d,i) {
if (d == -1) {
return "#000000";
} else {
return color(d);
}
});
}
function draw_legend() {
d3.select("#legend-info").selectAll("svg").remove();
var legend = Legend(d3.scaleLinear([0, max_concentration], ["#d0d0ec", "#7f0000"]), {
title: "Concentration",
tickFormat: ".2f",
tickValues: [0, max_concentration/4, max_concentration/2, max_concentration/4*3, max_concentration]
});
d3.select("#legend-info")
.node()
.appendChild(legend);
}
// tooltip events
var mouseover = function(event, d) {
tooltip.style("opacity", 1);
}
var mousemove = function(event, d) {
const tooltip_text = d3.select('.tooltip-text');
const tooltip_box = d3.select('.tooltip-box');
tooltip_text.text(`Value: ${d.value}`);
var box_width = tooltip_text.node().getComputedTextLength();
tooltip_box.attr("width", box_width+10);
const [x, y] = d3.pointer(event);
tooltip
.attr('transform', `translate(${x+svgMargin.left+50}, ${y+svgMargin.top})`);
}
var mouseleave = function(event, d) {
tooltip.style("opacity", 0)
}
// key events for changing grid iteration state
addEventListener("keydown", (event) => {
if (event.isComposing || event.keyCode === 81) {
if (state > 0) {
state -= 1;
max_concentration = calc_max_concentration(state);
create_n_update_scales();
draw_legend();
draw_n_update_grid(data[state]);
}
}
if (event.isComposing || event.keyCode === 69) {
if (state < iteration_count-1) {
state += 1;
max_concentration = calc_max_concentration(state);
create_n_update_scales();
draw_legend();
draw_n_update_grid(data[state]);
}
}
if (event.isComposing || event.keyCode === 65) {
if (state > 9) {
state -= 10;
} else {
state = 0;
}
max_concentration = calc_max_concentration(state);
create_n_update_scales();
draw_legend();
draw_n_update_grid(data[state]);
}
if (event.isComposing || event.keyCode === 68) {
if (state < iteration_count-10) {
state += 10;
} else {
state = iteration_count-1;
}
max_concentration = calc_max_concentration(state);
create_n_update_scales();
draw_legend();
draw_n_update_grid(data[state]);
}
});
// read data from file when selecting new file
async function create_data_struct(csv_text) {
var row_data = d3.dsvFormat(" ").parseRows(csv_text);
leftb = row_data[0];
rightb = row_data[1];
topb = row_data[2];
bottomb = row_data[3];
number_rows = row_data[0].length;
number_cols = row_data[2].length;
iteration_count = (row_data.length - 6) / (number_rows + 2);
console.log(iteration_count);
var iteration = []; // loop temp variable
var entry = {}; // loop temp variable
for (let i = 0; i < iteration_count; i++) {
iteration = []; // reset temp variable
for (let j = 0; j < number_rows; j++) {
for (let k = 0; k < number_cols; k++) {
entry = {}; // reset temp variable
entry.row = j;
entry.col = k;
entry.value = row_data[6 + i*(number_rows+2) + j][k];
iteration.push(entry);
}
}
data.push(iteration);
}
}
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", (evt) => {
if (!evt.target.files) { // no file, do nothing
return;
}
reader(evt.target.files[0], (err, res) => {
create_data_struct(res);
state = 0;
max_concentration = calc_max_concentration(state);
create_n_update_scales();
draw_legend();
draw_boundaries();
draw_n_update_grid(data[0]);
});
});
const svg = d3.create("svg")
.attr("width", width)
.attr("height", height);
var left_side = svg.append("g")
.attr("class", "boundary")
.attr("transform", `translate(0, ${svgMargin.top+boundary_gap})`);
var right_side = svg.append("g")
.attr("class", "boundary")
.attr("transform", `translate(${width-svgMargin.right}, ${svgMargin.top+boundary_gap})`);
var top_side = svg.append("g")
.attr("class", "boundary")
.attr("transform", `translate(${svgMargin.left+boundary_gap}, 0)`);
var bottom_side = svg.append("g")
.attr("class", "boundary")
.attr("transform", `translate(${svgMargin.left+boundary_gap}, ${height-svgMargin.bottom})`);
var grid = svg.append("g")
.attr("class", "grid")
.attr("transform", `translate(${svgMargin.left + boundary_gap}, ${svgMargin.top + boundary_gap})`);
// tooltip
var tooltip = svg.append("g")
.attr("class", "tooltip-area")
.style("opacity", 0);
var tooltip_box = tooltip.append("rect")
.attr("class", "tooltip-box")
.attr("height", 25)
.attr("transform", `translate(-5,-18)`)
.attr("fill", "white")
.attr("stroke", "black");
tooltip.append("text")
.attr("class", "tooltip-text")
.text("Value: ");
content.append(svg.node());
</script>