JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Improvements to mesh and double mesh #13

Merged
JorgeGonz merged 6 commits from feature/doubleMesh into development 2021-04-05 09:49:26 +02:00
Conversation 3 Commits 6 Files changed 52 +10451 -2017
13 changed files with 10451 additions and 2017 deletions
Download patch file Download diff file Expand all files Collapse all files
Showing only changes of commit cbcefb06c8 - Show all commits

Final implementation of a specific mesh for MCC, added a case for

Prev Next 
cylFlow that used two meshes and the User Manual has been updated.
Jorge Gonzalez 2021-04-05 09:45:57 +02:00
commit cbcefb06c8

10
doc/user-manual/bibliography.bib
View file

@ -8,6 +8,16 @@
pages = {3--67},
}
@article{higginson2020corrected,
title={A corrected method for Coulomb scattering in arbitrarily weighted particle-in-cell plasma simulations},
author={Higginson, Drew Pitney and Holod, Ihor and Link, Anthony},
journal={Journal of Computational Physics},
volume={413},
pages={109450},
year={2020},
publisher={Elsevier}
}
@Misc{gfortranURL,
author = {GNU Project},
title = {gfortran - the GNU Fortran compiler},

BIN
doc/user-manual/fpakc_UserManual.pdf
View file

Binary file not shown.

184
doc/user-manual/fpakc_UserManual.tex
View file

@ -21,6 +21,10 @@
Author = {Jorge Gonzalez}
}
}
% Allows breaking of URL in bibliography.
\setcounter{biburllcpenalty}{7000}
\setcounter{biburlucpenalty}{8000}
\makeglossaries
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -31,6 +35,8 @@
\newacronym{cpu}{CPU}{Central Processing Unit}
\newacronym{json}{JSON}{JavaScript Object Notation}
\newacronym{io}{I/O}{input/output}
\newacronym{mcc}{MCC}{Monte-Carlo Collisions}
\newacronym{cs}{CS}{Coulomb Scattering}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newglossaryentry{openmp}{name={OpenMP},description={Shared-memory parallelization}}
@ -60,7 +66,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Introduction}
\section{About fpakc}
The \Gls{fpakc} is a simulation tool that models species in plasmas (ions, electrons and neutrals) following the trajectories of macro-particles as they move and interact between them and the boundaries of the domain.
The \Gls{fpakc} is a simulation tool that models species in plasma (ions, electrons and neutrals) following the trajectories of macro-particles as they move and interact between them and the boundaries of the domain.
Particles properties are scattered into a finite element mesh in 1, 2 or three dimensions, with the possibility to choose different geometries.
The official repository can be found at: \url{https://gitlab.com/JorgeGonz/fpakc.git}.
The code is currently in very early steps of development and further improvements are expected very soon.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -78,11 +86,11 @@
This ease the process of fixing bugs and improving the codes by a large team of developers.
For more information, please refer to the \acrshort{fpakc} Coding Style document.
\item \acrshort{fpakc} requires to be ease to use.
Input files are required to be in a \textit{human} format, meaning that the different options can be easily understander without constant reference to the user guide.
\acrshort{fpakc} is aimed to be used in a wide range of applications and by a variety of scientist: from very established ones to newcomers to the field and students.
Input files are required to be in a \textit{human} format, meaning that the different options can be easily understood without constant reference to the user guide.
\acrshort{fpakc} is aimed to be used in a wide range of applications and by a variety of scientist: from very established ones to newcomers to the field and also students.
\end{enumerate}
These are foundation stones of the code and code development and should always be followed, at least for the releases in the official repository.
These are foundation stones of \acrshort{fpakc} and its development and should always be followed, at least for the releases in the official repository.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{How to collaborate}
@ -92,22 +100,32 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Operation Scheme}
\section{The Particle Method}
\acrshort{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals).
These macro-particles represent a large amount of real particles.
\Gls{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals).
These macro-particles could represent a large amount of real particles.
For now own, macro-particles will be referred as just particles by abusing of language.
In the evolution of these particles, external forces (as the electromagnetic field), interaction between particles (as collisions) and interaction with the boundaries of the domain are included.
During the initiation phase, the input and mesh file(s) are reading.
If an initial distribution for a species is specified in the input file, particles to match that distribution are loaded into the cells.
At each time step, particles are first pushed accounting for possible acceleration by external forces.
Then, the cell in which the particle ends up is located.
If a boundary is encountered, the interaction between the particle and the boundary is calculated.
Next, collisions for the particles in each cell are carried on.
This may include different collision processes for each particle.
Finally, the particles properties are scattered into the mesh nodes.
These properties are density, momentum and the stress tensor.
Non-dimensional units are used for this, but output files are converted into dimensional units.
If requested, the electromagnetic field is computed.
The general steps performed in each iteration are:
\begin{enumerate}
\item Firstly, new particles are introduced into the domain as specified in the input file.
\item Particles are then pushed accounting for possible acceleration by external forces.
During this process, if a particle changes cell it is found using the connectivity between elements.
If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall is computed.
A particle may abandon the computational domain and is no longer accounted for.
\item Next, collisions for the particles inside each cell are carried out.
This may include different collision processes for each particle.
Monte-Carlo collisions (elastic, ionization, charge-exchange$\ldots$) can be carried out in a specific mesh, to better adjust to the cell size required, similar to the mean-free path.
Although not yet implemented, Coulomb scattering will always be performed in the mesh used for scattering, which cell size should be in the order of the Debye length.
\item A new array containing all particles inside the numerical domain is obtained.
\item Finally, particle properties are scattered among the mesh nodes.
These properties are density, momentum and the stress tensor.
\item If requested, the electromagnetic field is computed.
\item If the number of iteration requires writing output files, it is done after all steps for the particles is completed.
\end{enumerate}
More in depth explanation of the different steps are given in the following sections.
\Gls{fpakc} has the capability to configure all the behavior of the simulation via the input file.
Parameters as injection, the kind of pusher used for each species, boundary conditions or collisions are user-input parameters and will be described in Chap.~\ref{ch:input_file}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Injection of new particles}
@ -122,43 +140,86 @@
Particles are pushed in the selected domain.
Velocity and position are updated according to the old particle values and the external forces.
All the push routines for the different geometries can be found in \lstinline|moduleSolver|.
The pushers included in \Gls{fpakc} are:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{3D Cartesian pusher}
Moving particles in a simple 3D Cartesian space.
\begin{itemize}
\item 3D Cartesian pusher.
Moving particles in a simple 3D Cartesian space.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{2D Cylindrical}
When a 2D cylindrical geometry is used ($z$, $r$), a Boris solver\cite{boris1970relativistic} is used to move particles accounting for the effect of the symmetry axis.
This pusher removes the issue with particles going to infinite velocity when $r \rightarrow 0$ by pushing the particles in Cartesian space and then converting it to $r-z$ geometry.
Velocity in the $\theta$ direction is updated for collision processes, although the dynamic in the angular direction is assumed as symmetric.
\item 2D Cylindrical.
When a 2D cylindrical geometry is used ($z$, $r$), a Boris solver\cite{boris1970relativistic} is used to move particles accounting for the effect of the symmetry axis.
This pusher removes the issue with particles going to infinite velocity when $r \rightarrow 0$ by pushing the particles in Cartesian space and then converting it to $r-z$ geometry.
Velocity in the $\theta$ direction is updated for collision processes, although the dynamic in the angular direction is assumed as symmetric.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{2D Cartesian pusher}
Moving particles in a simple 2D Cartesian space.
\item 2D Cartesian pusher.
Moving particles in a simple 2D Cartesian space.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{1D Radial pusher}
Same implementation as 2D cylindrical pusher but direction $z$ is ignored.
\item 1D Radial pusher.
Same implementation as 2D cylindrical pusher but direction $z$ is ignored.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{1D Cartesian pusher}
Moving particles in a simple 1D Cartesian space. Same implementation as in 2D Cartesian but $y$ direction is ignored.
\item 1D Cartesian pusher.
Moving particles in a simple 1D Cartesian space. Same implementation as in 2D Cartesian but $y$ direction is ignored.
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Find new cell}
Once the position and velocity of the particle are updated, the new cell that contains the particle is searched.
This is done by a neighbor search, starting from the previous cell containing the particle.
In the process of finding the new cell, it is possible that a particle encounters a boundary.
When the particle interacts with the boundary, the particle may continue its life in the simulation (reflected) or might be eliminated from it (absorbed).
When the particle interacts with the boundary, the particle may continue its life in the simulation or might be eliminated from it.
Once that the new cell is found or that the particle life has been terminated, the pushing is complete.
If a secondary mesh is used for the Monte-Carlo Collision method, the new cell in that mesh in which the particle reside is also found by the same method, although no interaction with the boundaries is accounted for this step.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Variable Weighting Scheme\label{sec:weightingScheme}}
One of the issues in particle simulations, specially for axial-symmetrical cases, is that due to the disparate volume of cells, specially close to the axis, the statistics in some cells is usually poor.
To try to fix that, the possibility to include a Variable Weighting Scheme in the simulations is available in \Gls{fpakc}.
This schemes detect when a particle change cells and modify its weight accordingly.
These schemes detect when a particle change cells and modify its weight accordingly.
To avoid particles having a larger weight than the rest, particle can be split in multiple particles if weight become too large.
The use of a Variable Weighting Scheme is defined by the user in the input file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Interaction between species}\label{ssec:collisions}
For each cell, interaction among the particles in it are carried out.
\Gls{fpakc} distinguish between two types of interactions: \acrfull{mcc} and \acrfull{cs}.
\acrshort{mcc} refers to the process in which two particles interact in short range.
These processes include, but are not limited to: elastic collisions, ionization/recombination, charge-exchange, excitation/de-excitation\ldots
A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collisions.
\acrshort{cs} refers to the large range interaction that a charged species suffer do to the charge of other particles.
The interactions between the different species is defined by the user.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{\acrlong{mcc}}
For each cell the maximum number of collisions between particle is computed.
For each collision, a random pair of particles is chosen.
A loop over all possible collisions for the pair of particles chosen is performed.
If a random number is above the probability of collision for that specific type, the collision take place.
If not, the next type for the particle pair is checked.
Below are described the type of collision process implemented in \acrshort{fpakc}:
\begin{itemize}
\item Elastic.
In this type of collision, particles exchange energy due to hard-sphere model.
Total energy is conserved.
Resulting velocity directions are chosen from Maxwellian distribution functions.
This interaction is useful for short-range collisions as neutral-neutral and charged-neutral elastic collisions.
\item Charge Exchange.
When an ion interacts with a neutral particle, an electron is exchanged between the two particles with no exchange of energy.
This is called a resonant charge-exchange.
\item Electron Impact Ionization.
When the relative energy between a neutral and an electron is above the ionization threshold, there is a probability that the neutral particle will become ionized.
This ionization emits and additional electron
\item Recombination.
When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
The photon emitted by this process is not modeled yet.
\end{itemize}
\subsection{\acrlong{cs}}
Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Reset of particle array}
@ -167,41 +228,6 @@
In this section, particles are assigned to the list of particles inside each individual cell.
Unfortunately, this is done right now without parallelisation and is very CPU consuming.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Interaction between species}\label{ssec:collisions}
For each cell, interaction among the particles in it are carried out.
The type of interaction between the different particles is defined by the user.
In general, the maximum number of interaction in a cell is computed.
For each collision, a pair of particles is selected.
A loop over all possible collisions for the pair of particles is performed.
If a random number generated is above the probability of collision for the type divided by the maximum one, the collision take place.
Collisions can change the velocity of the particles involved (elastic), create new particles (ionization-recombination) or change the type of particle (charge-exchange).
Below are described the type of collision process implemented between species.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elastic collision}
In this type of collision, particles exchange energy due to hard-sphere model.
Total energy is conserved.
Resulting velocity directions are chosen from Maxwellian distribution functions.
This interaction is useful for short-range collisions as neutral-neutral and charged-neutral elastic collisions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Charge Exchange}
When an ion interacts with a neutral particle, an electron is exchanged between the two particles with no exchange of energy.
This is called a resonant charge-exchange.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Electron Impact Ionization}
When the relative energy between a neutral and an electron is above the ionization threshold, there is a probability that the neutral particle will become ionized.
This ionization emits and additional electron
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Recombination}
When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
The photon emitted by this process is not modeled yet.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Scattering}
The properties of each particle are deposited in the nodes from the containing cell.
@ -620,9 +646,14 @@ make
\begin{itemize}
\item \textbf{folderCollisions}: Character.
Indicates the path to in which the cross section tables are allocated.
\item \textbf{meshCollisions}: Character.
Determines a specific mesh for \acrshort{mcc} processes.
The file needs to be located in the folder \textbf{output.folder}.
If this value is not present, the mesh defined in \textbf{geometry.meshFile} is used for \acrshort{mcc}.
The format of this mesh needs to be the same as the one defined in \textbf{geometry.meshType}.
\item \textbf{collisions}: Object.
Array.
Contains the different binary collisions.
Contains the different short range interactions (\acrshort{mcc}).
Multiple collision types can be defined for each pair of species.
Each object in the array is defined by:
\begin{itemize}
@ -665,21 +696,22 @@ make
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Example runs\label{ch:exampleRuns}}
\section{1D Cathode}
\section{1D Emissive Cathode (1D\_Cathode)}
Emission from a 1D cathond in both, cartesian and radial coordinates.
Both cases insert the same amount of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifiying some parameters in the input file.
The same mesh file (\lstinline|mesh.msh|) is used for both cases.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ALPHIE Grid system}
\section{ALPHIE Grid system (ALPHIE\_Grid)}
Two-dimensional axialsymmetry case to study the counterflow of electrons and Argon ions going trhough the ALPHIE grid system.
A \lstinline|mesh.geo| file is provided to easily modify the parameters of the grid system and generate a new mesh with \Gls{gmsh}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Flow around cylinder}
\section{Flow around cylinder (cylFlow)}
Simple case of neutral Argon flow around a cylinder in a 2D axialsymmetry geometry.
Elastic collisions between argon particles are included as default.
Elastic collisions between argon particles are included.
Two cases are presented here: one in which the same mesh (meshSingle.msh) for scattering and collisions is used (input.json) and a second one (inputDualMesh.json) in which a mesh is used for scattering (mesh.msh) and a second one is used only for collisions (meshColl.msh).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\printglossaries

4
runs/cylFlow/input.json
View file

@ -8,7 +8,7 @@
"geometry": {
"type": "2DCyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
"meshFile": "meshSingle.msh"
},
"species": [
{"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 5.0e8}
@ -57,7 +57,7 @@
},
"parallel": {
"OpenMP":{
"nThreads": 8
"nThreads": 24
}
}
}

64
runs/cylFlow/inputDualMesh.json Normal file
View file

@ -0,0 +1,64 @@
{
"output": {
"path": "./runs/cylFlow/",
"triggerOutput": 10,
"cpuTime": true,
"numColl": true
},
"geometry": {
"type": "2DCyl",
"meshType": "gmsh2",
"meshFile": "mesh.msh"
},
"species": [
{"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 5.0e8}
],
"boundary": [
{"name": "Injection", "physicalSurface": 1, "bTypes": [
{"type": "transparent"}
]},
{"name": "Chamber Walls", "physicalSurface": 2, "bTypes": [
{"type": "reflection"}
]},
{"name": "Exterior", "physicalSurface": 3, "bTypes": [
{"type": "transparent"}
]},
{"name": "Cylinder Walls", "physicalSurface": 4, "bTypes": [
{"type": "reflection"}
]},
{"name": "Axis", "physicalSurface": 5, "bTypes": [
{"type": "axis"}
]}
],
"inject": [
{"name": "Nozzle", "species": "Argon", "flow": 10.0, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1}
],
"reference": {
"density": 1.0e20,
"mass": 6.633e-26,
"temperature": 300.0,
"radius": 1.88e-10
},
"case": {
"tau": [5.0e-7],
"time": 1.0e-3,
"pusher": ["2DCylNeutral"],
"WeightingScheme": "Volume"
},
"interactions": {
"folderCollisions": "./data/collisions/",
"meshCollisions": "meshColl.msh",
"collisions": [
{"species_i": "Argon", "species_j": "Argon",
"cTypes": [
{"type": "elastic", "crossSection": "EL_Ar-Ar.dat"}
]}
]
},
"parallel": {
"OpenMP":{
"nThreads": 24
}
}
}

9
runs/cylFlow/mesh.geo
View file

@ -62,10 +62,13 @@ Physical Surface(4) = {4};
Physical Surface(5) = {5};
Transfinite Line {12, 2, 4, 6} = cyl_h/Lcell + 1 Using Progression 1;
Transfinite Line {1, 13, 10} = cyl_s/Lcell + 1 Using Progression 1;
Transfinite Line {11, 16, 15, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
Transfinite Line {10} = cyl_s/Lcell + 1 Using Progression 1;
Transfinite Line {1, 13} = cyl_s/Lcell + 1 Using Progression 0.95;
Transfinite Line {11, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
Transfinite Line {16,-15} = (dom_h - cyl_h)/Lcell + 1 Using Progression 0.95;
Transfinite Line {3, 9} = cyl_l/Lcell + 1 Using Progression 1;
Transfinite Line {5, 14, 8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
Transfinite Line {8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
Transfinite Line {-5, -14} = (dom_l - cyl_e)/Lcell + 1 Using Progression 0.95;
Transfinite Surface{1};
Recombine Surface {1};

3750
runs/cylFlow/mesh.msh
View file

File diff suppressed because it is too large Load diff

73
runs/cylFlow/meshColl.geo Normal file
View file

@ -0,0 +1,73 @@
cl__1 = 1;
cyl_h = 0.005;
cyl_l = 0.02;
cyl_s = 0.03;
cyl_e = cyl_s + cyl_l;
dom_h = 0.03;
dom_l = 0.07;
Lcell = 0.001;
Point(1) = {0, 0, 0, cl__1};
Point(2) = {cyl_s, 0, 0, cl__1};
Point(3) = {cyl_s, cyl_h, 0, cl__1};
Point(4) = {cyl_e, cyl_h, 0, cl__1};
Point(5) = {cyl_e, 0, 0, cl__1};
Point(6) = {dom_l, 0, 0, cl__1};
Point(7) = {dom_l, cyl_h, 0, cl__1};
Point(8) = {dom_l, dom_h, 0, cl__1};
Point(9) = {cyl_e, dom_h, 0, cl__1};
Point(10) = {cyl_s, dom_h, 0, cl__1};
Point(11) = {0, dom_h, 0, cl__1};
Point(12) = {0, cyl_h, 0, cl__1};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9};
Line(9) = {9, 10};
Line(10) = {10, 11};
Line(11) = {11, 12};
Line(12) = {12, 1};
Line(13) = {12, 3};
Line(14) = {4, 7};
Line(15) = {4, 9};
Line(16) = {10, 3};
Line Loop(1) = {1, 2, -13, 12};
Plane Surface(1) = {1};
Line Loop(2) = {13, -16, 10, 11};
Plane Surface(2) = {2};
Line Loop(3) = {3, 15, 9, 16};
Plane Surface(3) = {3};
Line Loop(4) = {5, 6, -14, 4};
Plane Surface(4) = {4};
Line Loop(5) = {14, 7, 8, -15};
Plane Surface(5) = {5};
Physical Surface(1) = {1};
Physical Surface(2) = {2};
Physical Surface(3) = {3};
Physical Surface(4) = {4};
Physical Surface(5) = {5};
Transfinite Line {12, 2, 4, 6} = cyl_h/Lcell + 1 Using Progression 1;
Transfinite Line {1, 13, 10} = cyl_s/Lcell + 1 Using Progression 1;
Transfinite Line {11, 16, 15, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
Transfinite Line {3, 9} = cyl_l/Lcell + 1 Using Progression 1;
Transfinite Line {5, 14, 8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
Transfinite Surface{1};
Recombine Surface {1};
Transfinite Surface{2};
Recombine Surface {2};
Transfinite Surface{3};
Recombine Surface {3};
Transfinite Surface{4};
Recombine Surface {4};
Transfinite Surface{5};
Recombine Surface {5};

3974
runs/cylFlow/meshColl.msh Normal file
View file

File diff suppressed because it is too large Load diff

79
runs/cylFlow/meshSingle.geo Normal file
View file

@ -0,0 +1,79 @@
cl__1 = 1;
cyl_h = 0.005;
cyl_l = 0.02;
cyl_s = 0.03;
cyl_e = cyl_s + cyl_l;
dom_h = 0.03;
dom_l = 0.07;
Lcell = 0.001;
Point(1) = {0, 0, 0, cl__1};
Point(2) = {cyl_s, 0, 0, cl__1};
Point(3) = {cyl_s, cyl_h, 0, cl__1};
Point(4) = {cyl_e, cyl_h, 0, cl__1};
Point(5) = {cyl_e, 0, 0, cl__1};
Point(6) = {dom_l, 0, 0, cl__1};
Point(7) = {dom_l, cyl_h, 0, cl__1};
Point(8) = {dom_l, dom_h, 0, cl__1};
Point(9) = {cyl_e, dom_h, 0, cl__1};
Point(10) = {cyl_s, dom_h, 0, cl__1};
Point(11) = {0, dom_h, 0, cl__1};
Point(12) = {0, cyl_h, 0, cl__1};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9};
Line(9) = {9, 10};
Line(10) = {10, 11};
Line(11) = {11, 12};
Line(12) = {12, 1};
Line(13) = {12, 3};
Line(14) = {4, 7};
Line(15) = {4, 9};
Line(16) = {10, 3};
Line Loop(1) = {1, 2, -13, 12};
Plane Surface(1) = {1};
Line Loop(2) = {13, -16, 10, 11};
Plane Surface(2) = {2};
Line Loop(3) = {3, 15, 9, 16};
Plane Surface(3) = {3};
Line Loop(4) = {5, 6, -14, 4};
Plane Surface(4) = {4};
Line Loop(5) = {14, 7, 8, -15};
Plane Surface(5) = {5};
Physical Line(1) = {12, 11};
Physical Line(2) = {10, 9, 8};
Physical Line(3) = {7, 6};
Physical Line(4) = {2, 3, 4};
Physical Line(5) = {1, 5};
Physical Surface(1) = {1};
Physical Surface(2) = {2};
Physical Surface(3) = {3};
Physical Surface(4) = {4};
Physical Surface(5) = {5};
Transfinite Line {12, 2, 4, 6} = cyl_h/Lcell + 1 Using Progression 1;
Transfinite Line {1, 13, 10} = cyl_s/Lcell + 1 Using Progression 1;
Transfinite Line {11, 16, 15, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
Transfinite Line {3, 9} = cyl_l/Lcell + 1 Using Progression 1;
Transfinite Line {5, 14, 8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
Transfinite Surface{1};
Recombine Surface {1};
Transfinite Surface{2};
Recombine Surface {2};
Transfinite Surface{3};
Recombine Surface {3};
Transfinite Surface{4};
Recombine Surface {4};
Transfinite Surface{5};
Recombine Surface {5};

4182
runs/cylFlow/meshSingle.msh Normal file
View file

File diff suppressed because it is too large Load diff

53
src/modules/mesh/moduleMesh.f90
View file

@ -400,6 +400,33 @@ MODULE moduleMesh
END INTERFACE
CONTAINS
!Constructs the global K matrix
SUBROUTINE constructGlobalK(self)
IMPLICIT NONE
CLASS(meshParticles), INTENT(inout):: self
INTEGER:: e
INTEGER, ALLOCATABLE:: n(:)
REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j
DO e = 1, self%numVols
n = self%vols(e)%obj%getNodes()
localK = self%vols(e)%obj%elemK()
nNodes = SIZE(n)
DO i = 1, nNodes
DO j = 1, nNodes
self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
END DO
END DO
END DO
END SUBROUTINE constructGlobalK
!Reset the output of node
PURE SUBROUTINE resetOutput(self)
USE moduleSpecies
@ -627,31 +654,11 @@ MODULE moduleMesh
END SUBROUTINE doCollisions
!Constructs the global K matrix
SUBROUTINE constructGlobalK(self)
IMPLICIT NONE
SUBROUTINE doCoulomb(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(inout):: self
INTEGER:: e
INTEGER, ALLOCATABLE:: n(:)
REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j
DO e = 1, self%numVols
n = self%vols(e)%obj%getNodes()
localK = self%vols(e)%obj%elemK()
nNodes = SIZE(n)
DO i = 1, nNodes
DO j = 1, nNodes
self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
END DO
END DO
END DO
END SUBROUTINE constructGlobalK
END SUBROUTINE doCoulomb
END MODULE moduleMesh

10
src/modules/moduleMath.f90
View file

@ -39,4 +39,14 @@ MODULE moduleMath
END FUNCTION normalize
!Norm 1 of vector
PURE FUNCTION norm1(a) RESULT(b)
IMPLICIT NONE
REAL(8), DIMENSION(:), INTENT(in):: a
REAL(8):: b
b = SUM(DABS(a))
END FUNCTION norm1
END MODULE moduleMath