Commit d18a2347 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

tidy up of input files. Some minor documentation


Former-commit-id: cfb30e69
parent 394af8a5
......@@ -6,13 +6,13 @@
class Checkerboard_SN : public ProblemBase
{
private:
Vector _scatteringXS; /*! @brief Vector of scattering crosssections */
Vector _scatteringXS; /*! @brief Vector of scattering crosssections */
Vector _totalXS; /*! @brief Vector of total crosssections */
Checkerboard_SN() = delete;
bool isAbsorption( const Vector& pos ) const; /*! @return True if pos is in absorption region, False otherwise */
bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
public:
Checkerboard_SN( Config* settings, Mesh* mesh );
......@@ -27,13 +27,13 @@ class Checkerboard_SN : public ProblemBase
class Checkerboard_PN : public ProblemBase
{
private:
Vector _scatteringXS; /*! @brief Vector of scattering crosssections */
Vector _totalXS; /*! @brief Vector of total crosssections */
Vector _scatteringXS; /*! @brief Vector of scattering crosssections len: numCells */
Vector _totalXS; /*! @brief Vector of total crosssections. len: numCells*/
Checkerboard_PN() = delete;
bool isAbsorption( const Vector& pos ) const; /*! @return True if pos is in absorption region, False otherwise */
bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
/**
* @brief Gets the global index for given order l of Legendre polynomials and given
......
......@@ -32,11 +32,11 @@ class Solver
double _dE; /*! @brief energy/time step size */
Vector _energies; // energy groups used in the simulation [keV]
std::vector<double> _density; // patient density, dim(_density) = _nCells
Vector _s; // stopping power, dim(_s) = _nTimeSteps
Vector _s; // stopping power, dim(_s) = _nTimeSteps (only for csdsolver)
std::vector<VectorVector> _Q; /*! @brief external source term */
VectorVector _sigmaS; /*! @brief scattering cross section for all energies */
VectorVector _sigmaT; /*! @brief total cross section for all energies */
VectorVector _sigmaS; /*! @brief scattering cross section for all energies. len: _nEnergies x numCells */
VectorVector _sigmaT; /*! @brief total cross section for all energies. len: _nEnergies x numCells*/
// quadrature related numbers
QuadratureBase* _quadrature; /*! @brief quadrature to create members below */
......
......@@ -57,7 +57,7 @@ inline void PrintVectorVector( const VectorVector vectorIn ) {
else
printf( "| %.2f ", vectorIn[idx_Outer][idx_Inner] );
}
std::cout << " |\n";
printf( " |\n" );
}
}
......
// settings Kersting
lc = 0.005;
width = 1.0;
length = 10.0;
start = -2.5;
//+
Point(1) = {start, width, 0, lc};
//+
Point(2) = {start+length, width, 0, lc};
//+
Point(3) = {start, 0, 0, lc};
//+
Point(4) = {start+length, 0, 0, lc};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 4};
//+
Line(3) = {4, 3};
//+
Line(4) = {3, 1};
//+
Curve Loop(1) = {4, 1, 2, 3};
//+
Plane Surface(1) = {1};
//+
Physical Curve("dirichlet") = {4, 2};
//+
Physical Curve("wall_low") = {1};
Physical Curve("wall_up") = {3};
//+
Transfinite Surface {1} = {1, 2, 4, 3};
//+
Transfinite Curve {4, 2} = 1 Using Progression 1;
//+
Transfinite Curve {1, 3} = 1000 Using Progression 1;
//+
Recombine Surface {1};
This diff is collapsed.
% ---- File specifications ----
OUTPUT_DIR = ../result
OUTPUT_FILE = checkerboard
LOG_DIR = ../result/logs
MESH_FILE = checkerboard.su2
% ---- Solver specifications ----
CFL_NUMBER = 0.5
TIME_FINAL = 3.2
PROBLEM = CHECKERBOARD
% ---- Boundary Conditions ----
BC_NEUMANN = ( void )
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
QUAD_ORDER = 10
% ----- Output ----
%
VOLUME_OUTPUT = (MINIMAL)
VOLUME_OUTPUT_FREQUENCY = 1
SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
SCREEN_OUTPUT_FREQUENCY = 1
HISTORY_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
HISTORY_OUTPUT_FREQUENCY = 1
import pygmsh as pg
import numpy as np
import itertools
import os
def add_block(x0,y0,length,char_length,geom):
coords = np.array([
[x0, y0, 0.0],
[x0+length, y0, 0.0],
[x0+length, y0+length, 0.0],
[x0, y0+length, 0.0]
])
return geom.add_polygon(coords, char_length)
char_length = 0.075
geom = pg.opencascade.Geometry()
domain = add_block(0, 0, 7, char_length, geom)
xpos = ypos = [1, 2, 3, 4, 5]
pattern = list(itertools.product(xpos, ypos))[::2]
pattern.pop(7)
boxes = [domain]
for pos in pattern:
boxes.append(add_block(pos[0], pos[1], 1, char_length, geom))
geom.boolean_fragments(boxes,[])
geom.add_physical(domain.lines, label="void")
mesh_code = geom.get_code()
with open("checkerboard.geo","w") as mesh_file:
mesh_file.write(mesh_code)
os.system('gmsh checkerboard.geo -2 -format su2 -save_all')
os.system('rm checkerboard.geo')
This diff is collapsed.
......@@ -12,7 +12,7 @@ OUTPUT_DIR = ../result
OUTPUT_FILE = example_csd1
LOG_DIR = ../result/logs
% Mesh File
MESH_FILE = 1DMesh.su2
MESH_FILE = meshes/1DMesh.su2
%
PROBLEM = WATERPHANTOM
%
......
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data Generator for %
% Neural Entropy Closure %
% Author <Steffen Schotthöfer> %
% Date 17.12.2020 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ---- Global settings ----
DATA_GENERATOR_MODE = YES
TRAINING_SET_SIZE = 500
MAX_VALUE_FIRST_MOMENT = 20
REALIZABLE_SET_EPSILON_U0 = 0.05
REALIZABLE_SET_EPSILON_U1 = 0.97
%
% ---- File specifications ----
%
% Output directory
OUTPUT_DIR = ../result
% Output file
OUTPUT_FILE = trainM0
% Log directory
LOG_DIR = ../result/logs
%
% --- Spherical Basis ----
%
% Choice of basis functions
SPHERICAL_BASIS = SPHERICAL_MONOMIALS
%
% Maximal Moment degree
MAX_MOMENT_SOLVER = 1
%
% --- Entropy settings ----
ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
ENTROPY_OPTIMIZER = NEWTON
%
% ----- Newton Solver Specifications ----
%
NEWTON_FAST_MODE = NO
NEWTON_ITER = 1000000
NEWTON_EPSILON = 0.01
NEWTON_STEP_SIZE = 0.7
NEWTON_LINE_SEARCH_ITER = 1000
%
% ----- Quadrature ----
%
% Quadrature Rule
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
% Quadrature Order
QUAD_ORDER = 8
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linesource Benchmarking File MN %
% Author <Steffen Schotthöfer> %
% Date 01.07.2020 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% ---- File specifications ----
%
% Output directory
OUTPUT_DIR = ../result
% Output file
OUTPUT_FILE = exampleMN_mono
% Log directory
LOG_DIR = ../result/logs
% Mesh File
MESH_FILE = linesource.su2
%MESH_FILE = linesource_debug.su2
%
% ---- Problem description ---
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 1
%
% ---- Solver specifications ----
%
% Choice of basis functions
SPHERICAL_BASIS = SPHERICAL_MONOMIALS
% Maximal Moment degree
MAX_MOMENT_SOLVER = 1
% Solver type
SOLVER = MN_SOLVER
% CFL number
CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.5
% Reconstruction order (spatial flux)
RECONS_ORDER = 1
%
% ---- Entropy settings ----
%
ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
ENTROPY_OPTIMIZER = NEWTON
NEURAL_MODEL = 4
%
% ----- Newton Solver Specifications ----
%
NEWTON_FAST_MODE = NO
NEWTON_ITER = 1000
NEWTON_EPSILON = 0.01
NEWTON_STEP_SIZE = 0.7
NEWTON_LINE_SEARCH_ITER = 1000
%
% ---- Boundary Conditions ----
%
% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
% Dirichlet Boundary
BC_DIRICHLET = ( void )
%
% ----- Quadrature ----
%
% Quadrature Rule
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
% Quadrature Order
QUAD_ORDER = 8
%
% ----- Output ----
%
VOLUME_OUTPUT = (MINIMAL, MOMENTS, DUAL_MOMENTS)
VOLUME_OUTPUT_FREQUENCY = 1
SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
SCREEN_OUTPUT_FREQUENCY = 1
HISTORY_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
HISTORY_OUTPUT_FREQUENCY = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linesource Benchmarking File PN %
% Author <Steffen Schotthöfer> %
% Date 01.07.2020 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% ---- File specifications ----
%
% Output directory
OUTPUT_DIR = ../result
% Output file
OUTPUT_FILE = examplePN
% Log directory
LOG_DIR = ../result/logs
% Mesh File
MESH_FILE = linesource.su2
%MESH_FILE = linesource_debug.su2
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 1
%
%
% ---- Solver specifications ----
%
% Solver type
SOLVER = PN_SOLVER
% CFL number
CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.3
% Maximal Moment degree
MAX_MOMENT_SOLVER = 2
% Reconstruction order
RECONS_ORDER = 2
% ---- Boundary Conditions ----
% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
% Dirichlet Boundary
BC_DIRICHLET = ( void )
%
% ----- Output ----
%
VOLUME_OUTPUT = (MINIMAL)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linesource Benchmarking File SN %
% Author <Steffen Schotthöfer> %
% Date 01.10.2020 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% ---- File specifications ----
%
% Output directory
OUTPUT_DIR = ../result
% Output file
OUTPUT_FILE = example_SN
% Log directory
LOG_DIR = ../result/logs
% Mesh File
MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
%
% ---- Solver specifications ----
%
% Solver type
SOLVER = SN_SOLVER
% CFL number
CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.5
% Reconstruction order
RECONS_ORDER = 2
%
% ---- Boundary Conditions ----
% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
% Dirichlet Boundary
BC_DIRICHLET = ( void )
%
%% Quadrature Specifications
% Quadrature Type
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
% Quadrature Order
QUAD_ORDER = 8
%
% ----- Output ----
%
VOLUME_OUTPUT = (ANALYTIC, MINIMAL)
VOLUME_OUTPUT_FREQUENCY = 1
......@@ -4,7 +4,7 @@
OUTPUT_DIR = ../result
OUTPUT_FILE = example_csdSN_1d_NOScatter
LOG_DIR = ../result/logs
MESH_FILE = 1DMesh.su2
MESH_FILE = meshes/1DMesh.su2
PROBLEM = WATERPHANTOM
SOLVER = CSD_SN_SOLVER
......
import pygmsh as pg
import numpy as np
import itertools
import os
def add_block(x0,y0,lengthX, lengthY,char_length,geom):
coords = np.array([
[x0, y0, 0.0],
[x0+lengthX, y0, 0.0],
[x0+lengthX, y0+lengthY, 0.0],
[x0, y0+lengthY, 0.0]
])
return geom.add_polygon(coords, char_length)
char_length = 0.015
geom = pg.opencascade.Geometry()
domain = add_block(-1, -1, 2,2, char_length, geom)
geom.add_raw_code('psource = newp;\nPoint(psource) = {0.0, 0.0, 0.0, '+str(char_length)+'};\nPoint{psource} In Surface{'+domain.id+'};')
geom.add_physical(domain.lines, label="void")
mesh_code = geom.get_code()
with open("linesource.geo","w") as mesh_file:
mesh_file.write(mesh_code)
os.system('gmsh linesource.geo -2 -format su2 -save_all')
os.system('rm linesource.geo')
This diff is collapsed.
......@@ -95,7 +95,7 @@ std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*en
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( isSource( cellMids[j] ) ) Q[j] = 1.0 / std::sqrt( 4 * M_PI ); // isotropic source
if( isSource( cellMids[j] ) ) Q[j] = 10.0 / std::sqrt( 4 * M_PI ); // isotropic source
}
return std::vector<VectorVector>( 1u, Q );
}
......
......@@ -179,9 +179,6 @@ void MNSolver::FluxUpdate() {
// Dirichlet Boundaries stayd
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
_solNew[idx_cell] = ConstructFlux( idx_cell );
// if( norm( _solNew[idx_cell] ) > 0.001 ) {
// std::cout << _solNew[idx_cell] << "\n";
//}
}
}
......@@ -193,10 +190,9 @@ void MNSolver::FVMUpdate( unsigned idx_energy ) {
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// Flux update
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_sys] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys]; /* cell averaged flux */
-_dE* _sol[idx_cell][idx_sys] * ( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
}
// Source Term
_solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
......
......@@ -10,6 +10,7 @@
#include "solvers/csdsolvertrafofp.h"
#include "solvers/csdsolvertrafofp2d.h"
#include "solvers/csdsolvertrafofpsh2d.h"
#include "toolboxes/textprocessingtoolbox.h"
#include "solvers/mnsolver.h"
#include "solvers/pnsolver.h"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment