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

First implementation of probing method

Browse source 
The code nows offer the possibility to obtain the distribution function
for a specific species in a 3D velocity grid at a determined position.

This is a simple method that just scatter the particles in one cell into
the velocity grid.
This commit is contained in:
Jorge Gonzalez 2021-06-29 10:37:39 +02:00
commit 56967dd6c7
7 changed files with 413 additions and 86 deletions
Download patch file Download diff file Expand all files Collapse all files

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

Binary file not shown.

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

@ -215,11 +215,12 @@
\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.
The photon emitted by this process is not modelled yet.
\end{itemize}
\subsection{\acrlong{cs}}
Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \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}
@ -227,6 +228,20 @@
The new array containing only the particles inside the domain will be the one used in the next steps.
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{Probing}\label{sec:probing}
As default, \gls{fpakc} outputs information of macroscopic quantities (density, velocity, temperature\ldots) in the finite element mesh.
However, a lot of information can be extracted from the particle distribution function.
Thus, a probing method is provided to extract the distribution function in a specific position.
The particles inside a cell in which the input position is located are distributed into a 3D velocity grid.
The user can decide the grid width and the number of points in each direction.
The distribution function will be calculated and wrote with a time step decided by the user.
If a particle velocity resides outside of the velocity grid (in any direction), it wont be added to the tally of the distribution function.
Due to the limitation of only taking into account particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results.
A more advance method taking into account distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Scattering}
@ -359,6 +374,30 @@ make
Trigger between writings is the same as in \textbf{triggerOutput}.
\item \textbf{EMField}: Logical.
Determines if the electromagnetic field is printed.
\item \textbf{probes}: Array of objects.
Defines the probes employed for obtaining the distribution function at specific positions.
See Sec.~\ref{sec:probing} for more information.
The object is structured as follows:
\begin{itemize}
\item \textbf{species}: Character.
Species name as defined in \textbf{species} array.
\item \textbf{position}: Real.
Array of dimension $3$.
Units in $\unit{m}$.
Indicates the position of the probing.
\item \textbf{timeStep}: Real
Units in $\unit{s}$.
Optional.
Time step for output of the distribution function.
If none is provided, the minimum time step of the case is used.
\item \textbf{velocity\_1, velocity\_2, velocity\_3}: Real.
Array of dimension $2$.
Velocity range (minimum-maximum) in which the distribution function will be interpolated.
The subscripts 1, 2, 3 indicate the three directions of the case.
\item \textbf{points}: Integer.
Array of dimension $3$.
Number of points in each direction.
\end{itemize}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -439,72 +478,72 @@ make
\begin{itemize}
\item \textbf{name}: Character.
Name of the boundary.
\item \textbf{type}: Character.
Type of boundary.
Accepted values are:
\begin{itemize}
\item \textbf{reflection}: Elastic reflection of particles.
\item \textbf{absorption}: Particle is eliminated from the domain.
The particle is first moved into the edge and its properties are scattered among the edge nodes.
\item \textbf{transparent}: Particle abandon the numerical domain.
\item \textbf{wallTemperature}: Reflective wall with cosntant temperature that exchange heat with particles.
Required parameters are:
\begin{itemize}
\item \textbf{temperature}: Real.
Units of $\unit{K}$.
Temperature wall.
\item \textbf{specificHeat}: Real.
Units of $\unit{J kg^{-1} K^{-1}}$.
Specific heat capacity of the material.
\end{itemize}
\item \textbf{ionization}: Per each particle crossing the surface with this type of boundary, a number of ionization events are calculated.
A pair of ion-electron is generated for each ionization event taking as a reference a neutral background.
Secondary electron is taken as same type as incident particle.
The available input is:
\begin{itemize}
\item \textbf{neutral}: Object.
Information about the neutral background.
Required parameters are:
\begin{itemize}
\item \textbf{ion}: Character.
Species name of the ion generated as defined in object \textbf{species}.
Required parameter.
\item \textbf{mass}: Real.
Units in $\unit{kg}$.
Mass of neutral species.
If missing, the mass of the ion is ussed
\item \textbf{density}: Real.
Units in $\unit{m^{-3}}$.
Density of neutral background.
Required parameter.
\item \textbf{velocity}: Real.
Units in $\unit{m^{-3}}$.
Array of dimension $3$.
Mean velocity of neutral background.
Required parameter.
\item \textbf{temperature}: Real.
Units in $\unit{K}$.
Temperature of neutral background.
Required parameter.
\end{itemize}
\item \textbf{effectiveTime}: Real.
Units in $\unit{s}$.
As the particle is no longer simulated once it crossed the boundary, this time represent the effective time in which the particle produces ionization processes in the neutral background.
Required parameter.
\item \textbf{energyThreashold}: Real.
Units in $\unit{eV}$.
Ionization energy threshold for the simulated process.
Required parameter.
\item \textbf{crossSection}: Character.
Complete path to the cross section data for the ionization process.
\end{itemize}
\item \textbf{axis}: Identifies the symmetry axis for 2D cylindrical simulations.
If for some reason a particle interact with this axis, it is reflected.
\end{itemize}
\item \textbf{physicalSurface}: Integer.
Identification of the edge in the mesh file.
Identification of the surface in the mesh file.
\item \textbf{bType}: Array of objects of dimension 'number of species'.
Per each species defined in the case, a boundary \textbf{type} needs to be provided.
Accepted values for \textbf{type} are:
\begin{itemize}
\item \textbf{reflection}: Elastic reflection of particles.
\item \textbf{absorption}: Particle is eliminated from the domain.
The particle is first moved into the edge and its properties are scattered among the edge nodes.
\item \textbf{transparent}: Particle abandon the numerical domain.
\item \textbf{wallTemperature}: Reflective wall with cosntant temperature that exchange heat with particles.
Required parameters are:
\begin{itemize}
\item \textbf{temperature}: Real.
Units of $\unit{K}$.
Temperature wall.
\item \textbf{specificHeat}: Real.
Units of $\unit{J kg^{-1} K^{-1}}$.
Specific heat capacity of the material.
\end{itemize}
\item \textbf{ionization}: Per each particle crossing the surface with this type of boundary, a number of ionization events are calculated.
A pair of ion-electron is generated for each ionization event taking as a reference a neutral background.
Secondary electron is taken as same type as incident particle.
The available input is:
\begin{itemize}
\item \textbf{neutral}: Object.
Information about the neutral background.
Required parameters are:
\begin{itemize}
\item \textbf{ion}: Character.
Species name of the ion generated as defined in object \textbf{species}.
Required parameter.
\item \textbf{mass}: Real.
Units in $\unit{kg}$.
Mass of neutral species.
If missing, the mass of the ion is ussed
\item \textbf{density}: Real.
Units in $\unit{m^{-3}}$.
Density of neutral background.
Required parameter.
\item \textbf{velocity}: Real.
Units in $\unit{m^{-3}}$.
Array of dimension $3$.
Mean velocity of neutral background.
Required parameter.
\item \textbf{temperature}: Real.
Units in $\unit{K}$.
Temperature of neutral background.
Required parameter.
\end{itemize}
\item \textbf{effectiveTime}: Real.
Units in $\unit{s}$.
As the particle is no longer simulated once it crossed the boundary, this time represent the effective time in which the particle produces ionization processes in the neutral background.
Required parameter.
\item \textbf{energyThreashold}: Real.
Units in $\unit{eV}$.
Ionization energy threshold for the simulated process.
Required parameter.
\item \textbf{crossSection}: Character.
Complete path to the cross section data for the ionization process.
\end{itemize}
\item \textbf{axis}: Identifies the symmetry axis for 2D cylindrical simulations.
If for some reason a particle interact with this axis, it is reflected.
\end{itemize}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -631,8 +670,7 @@ make
\begin{itemize}
\item \textbf{Electrostatic}: Solves the Poison equation to obtain the self-consistent electrostatic potential.
\end{itemize}
\item \textbf{initial}: Object.
Array.
\item \textbf{initial}: Array of objects.
Determines initial values for the species.
Required values are:
\begin{itemize}
@ -661,8 +699,7 @@ make
Units of $\unit{s}$.
Time step to calculate MCC.
If no time is provided, the minimum time step is used.
\item \textbf{collisions}: Object.
Array.
\item \textbf{collisions}: Array of objects.
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:
@ -670,8 +707,7 @@ make
\item \textbf{species\_i}, \textbf{species\_j}: Character.
Define the two species involved in the collision processes.
Order is indiferent.
\item \textbf{cTypes}: Object.
Array.
\item \textbf{cTypes}: Array of objects.
Defines all the collisions between \textbf{species\_i} and \textbf{species\_j}.
Required values are:
\begin{itemize}
@ -705,29 +741,29 @@ make
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Example runs\label{ch:exampleRuns}}
\chapter{Example runs}\label{ch:exampleRuns}
This chapter presents a description of the different example files distributed with \acrshort{fpakc}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{1D Emissive Cathode (1D\_Cathode)}
Emission from a 1D cathond in both, cartesian and radial coordinates.
Emission from a 1D cathode 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.
This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file.
The same mesh file (\lstinline|mesh.msh|) is used for both cases.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{0D \ce{Ar}-\ce{Ar+} Elastic Collision (0D\_Argon)}
Test case to check the 0D geometry that includes the elastic collision between \ce{Ar} and \ce{Ar+}.
Initial states are readed from the Argon\_Initial.dat and Argon+\_Initial.dat
Initial states are read from the Argon\_Initial.dat and Argon+\_Initial.dat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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.
Two-dimensional axial-symmetry case to study the counterflow of electrons and Argon ions going through 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 (cylFlow)}
Simple case of neutral Argon flow around a cylinder in a 2D axialsymmetry geometry.
Simple case of neutral Argon flow around a cylinder in a 2D axial-symmetry geometry.
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).

10
src/fpakc.f90
View file

@ -1,12 +1,13 @@
! FPAKC main program
PROGRAM fpakc
USE moduleCompTime
USE moduleCaseParam
USE moduleInput
USE moduleErrors
USE moduleInject
USE moduleSolver
USE moduleMesh
USE moduleCompTime
USE moduleCaseParam
USE moduleProbe
USE moduleErrors
USE OMP_LIB
IMPLICIT NONE
@ -84,6 +85,9 @@ PROGRAM fpakc
!$OMP SINGLE
tCoul = omp_get_wTime() - tCoul
!Do probing
CALL doProbes(t)
!Reset particles
tReset = omp_get_wtime()
!$OMP END SINGLE

1
src/makefile
View file

@ -4,6 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
$(OBJDIR)/moduleProbe.o \
$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
$(OBJDIR)/moduleMesh3DCart.o \

5
src/modules/makefile
View file

@ -2,7 +2,8 @@
OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
moduleOutput.o moduleInput.o moduleSolver.o \
moduleCollisions.o moduleTable.o moduleParallel.o \
moduleEM.o moduleRandom.o moduleMath.o
moduleEM.o moduleRandom.o moduleMath.o \
moduleProbe.o
all: $(OBJS) mesh.o
@ -12,7 +13,7 @@ mesh.o: moduleCollisions.o moduleBoundary.o moduleMath.o
moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleInput.f90
moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleProbe.o moduleInput.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90

48
src/modules/moduleInput.f90
View file

@ -60,6 +60,11 @@ MODULE moduleInput
CALL readGeometry(config)
CALL checkStatus(config, "readGeometry")
!Read probes
CALL verboseError('Reading Probes...')
CALL readProbes(config)
CALL checkStatus(config, 'readProbes')
!Read initial state for species
CALL verboseError('Reading Initial state...')
CALL readInitial(config)
@ -903,6 +908,47 @@ MODULE moduleInput
END SUBROUTINE readGeometry
SUBROUTINE readProbes(config)
USE moduleProbe
USE moduleErrors
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(2):: istring
INTEGER:: i
CHARACTER(:), ALLOCATABLE:: speciesName
REAL(8), ALLOCATABLE, DIMENSION(:):: r
REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
INTEGER, ALLOCATABLE, DIMENSION(:):: points
REAL(8):: timeStep
CALL config%info('output.probes', found, n_children = nProbes)
IF (found) ALLOCATE(probe(1:nProbes))
DO i = 1, nProbes
WRITE(iString, '(I2)') i
object = 'output.probes(' // trim(istring) // ')'
CALL config%get(object // '.species', speciesName, found)
CALL config%get(object // '.position', r, found)
CALL config%get(object // '.velocity_1', v1, found)
CALL config%get(object // '.velocity_2', v2, found)
CALL config%get(object // '.velocity_3', v3, found)
CALL config%get(object // '.points', points, found)
CALL config%get(object // '.timeStep', timeStep, found)
IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
END DO
END SUBROUTINE readProbes
SUBROUTINE readEMBoundary(config)
USE moduleMesh
USE moduleOutput
@ -926,7 +972,7 @@ MODULE moduleInput
IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
DO i = 1, nBoundaryEM
WRITE(istring, '(i2)') i
WRITE(istring, '(I2)') i
object = 'boundaryEM(' // trim(istring) // ')'
CALL config%get(object // '.type', boundEM(i)%typeEM, found)

239
src/modules/moduleProbe.f90 Normal file
View file

@ -0,0 +1,239 @@
MODULE moduleProbe
USE moduleMesh
USE moduleSpecies
TYPE:: probeDistFunc
INTEGER:: id
REAL(8), DIMENSION(1:3):: r
CLASS(meshVol), POINTER:: cell
CLASS(speciesGeneric), POINTER:: species
REAL(8), ALLOCATABLE, DIMENSION(:):: vi, vj, vk
REAL(8), ALLOCATABLE, DIMENSION(:,:,:):: f
INTEGER, DIMENSION(1:3):: nv
REAL(8), DIMENSION(1:3):: vrange
REAL(8):: dvInv
INTEGER:: every
CONTAINS
PROCEDURE, PASS:: init
PROCEDURE, PASS:: findLowerIndex
PROCEDURE, PASS:: calculate
PROCEDURE, PASS:: output
END TYPE probeDistFunc
INTEGER:: nProbes = 0
TYPE(probeDistFunc), ALLOCATABLE:: probe(:)
CONTAINS
!Functions for probeDistFunc type
SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep)
USE moduleCaseParam
USE moduleRefParam
USE moduleSpecies
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(probeDistFunc), INTENT(out):: self
INTEGER, INTENT(in):: id
CHARACTER(:), ALLOCATABLE, INTENT(in):: speciesName
REAL(8), INTENT(in):: r(1:3)
REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
INTEGER, INTENT(in):: points(1:3)
REAL(8), INTENT(in):: timeStep
INTEGER:: sp, e, i
REAL(8):: dv(1:3)
!Assign id
self%id = id
!Get species
sp = speciesName2Index(speciesName)
self%species => species(sp)%obj
!Find cell
self%r = r/L_ref
e = findCellBrute(mesh, self%r)
IF (e == 0) CALL criticalError("No cell found for position in probe", 'init')
self%cell => mesh%vols(e)%obj
!Allocates velocity grid
self%nv = points
ALLOCATE(self%vi(1:self%nv(1)))
ALLOCATE(self%vj(1:self%nv(2)))
ALLOCATE(self%vk(1:self%nv(3)))
!Creates grid
dv(1) = (v1(2) - v1(1))/REAL(self%nv(1) - 1)/v_ref
dv(2) = (v2(2) - v2(1))/REAL(self%nv(2) - 1)/v_ref
dv(3) = (v3(2) - v3(1))/REAL(self%nv(3) - 1)/v_ref
self%dvInv = 1.D0/(dv(1)*dv(2)*dv(3))
DO i = 1, self%nv(1)
self%vi(i) = dv(1)*REAL(i - 1) + v1(1)/v_ref
END DO
DO i = 1, self%nv(2)
self%vj(i) = dv(2)*REAL(i - 1) + v2(1)/v_ref
END DO
DO i = 1, self%nv(3)
self%vk(i) = dv(3)*REAL(i - 1) + v3(1)/v_ref
END DO
self%vrange(1) = self%vi(self%nv(1)) - self%vi(1)
self%vrange(2) = self%vj(self%nv(2)) - self%vj(1)
self%vrange(3) = self%vk(self%nv(3)) - self%vk(1)
!Allocates distribution function
ALLOCATE(self%f(1:self%nv(1), &
1:self%nv(2), &
1:self%nv(3)))
!Number of iterations between output
self%every = NINT(timeStep/ tauMin / ti_ref)
END SUBROUTINE init
SUBROUTINE findLowerIndex(self, vp, i, j, k, inside)
IMPLICIT NONE
CLASS(probeDistFunc), INTENT(in):: self
REAL(8), INTENT(in):: vp(1:3)
INTEGER, INTENT(out):: i, j, k
LOGICAL, INTENT(out):: inside
i = FLOOR((vp(1) - self%vi(1))/self%vrange(1)*(REAL(self%nv(1) - 1)) + 1.D0)
IF (i >= self%nv(1) .OR. i < 1) inside = .FALSE.
j = FLOOR((vp(2) - self%vj(1))/self%vrange(2)*(REAL(self%nv(2) - 1)) + 1.D0)
IF (j >= self%nv(2) .OR. j < 1) inside = .FALSE.
k = FLOOR((vp(3) - self%vk(1))/self%vrange(3)*(REAL(self%nv(3) - 1)) + 1.D0)
IF (k >= self%nv(3) .OR. k < 1) inside = .FALSE.
END SUBROUTINE findLowerIndex
SUBROUTINE calculate(self)
USE moduleSpecies
USE moduleList
IMPLICIT NONE
CLASS(probeDistFunc), INTENT(inout):: self
TYPE(particle), POINTER:: part
TYPE(lNode), POINTER:: node
INTEGER:: i, j, k
LOGICAL:: inside
REAL(8):: fi, fi1
REAL(8):: fj, fj1
REAL(8):: fk, fk1
!Reset distribution function
self%f = 0.D0
!Loop over particles in cell
node => self%cell%listPart_in%head
DO WHILE(ASSOCIATED(node))
!Selects particle on list
part => node%part
!If particle is of desired type, include in the distribution function
IF (part%species%n == self%species%n) THEN
!find lower index for all dimensions
CALL self%findLowerIndex(part%v, i, j, k, inside)
!If particle is inside the velocity grid, add it to the distribution function
IF (inside) THEN
!Calculate weights
fi = self%vi(i+1) - part%v(1)
fi1 = part%v(1) - self%vi(i)
fj = self%vj(j+1) - part%v(2)
fj1 = part%v(2) - self%vj(j)
fk = self%vk(k+1) - part%v(3)
fk1 = part%v(3) - self%vk(k)
!Assign particle weight to distribution function
self%f(i , j , k ) = fi * fj * fk * part%weight
self%f(i+1, j , k ) = fi1 * fj * fk * part%weight
self%f(i , j+1, k ) = fi * fj1 * fk * part%weight
self%f(i+1, j+1, k ) = fi1 * fj1 * fk * part%weight
self%f(i , j , k+1) = fi * fj * fk1 * part%weight
self%f(i+1, j , k+1) = fi1 * fj * fk1 * part%weight
self%f(i , j+1, k+1) = fi * fj1 * fk1 * part%weight
self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * part%weight
END IF
END IF
!Move to next particle in the list
node => node%next
END DO
!Divide by the velocity cube volume
self%f = self%f * self%dvInv
END SUBROUTINE calculate
SUBROUTINE output(self, t)
USE moduleOutput
USE moduleRefParam
IMPLICIT NONE
CLASS(probeDistFunc), INTENT(in):: self
INTEGER, INTENT(in):: t
CHARACTER (LEN=iterationDigits):: tstring
CHARACTER (LEN=3):: pstring
CHARACTER(:), ALLOCATABLE:: filename
INTEGER:: i, j, k
WRITE(tstring, iterationFormat) t
WRITE(pstring, "(I3.3)") self%id
fileName='OUTPUT_' // tstring// '_f_' // pstring // '.dat'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (10, file = path // folder // '/' // fileName)
WRITE(10, "(A1, 1X, A)") "# ", self%species%name
WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(t)*tauMin*ti_ref, " s"
WRITE(10, "(A1, A19, 3(A20))") "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
DO i = 1, self%nv(1)
DO j = 1, self%nv(2)
DO k = 1, self%nv(3)
WRITE(10, "(4(ES20.6E3))") self%vi(i)*v_ref, &
self%vj(j)*v_ref, &
self%vk(k)*v_ref, &
self%f(i, j, k)
END DO
WRITE(10, *)
END DO
END DO
CLOSE(10)
END SUBROUTINE output
SUBROUTINE doProbes(t)
IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: i
DO i = 1, SIZE(probe)
IF (MOD(t, probe(i)%every) == 0 .OR. t == tFinal) THEN
CALL probe(i)%calculate()
CALL probe(i)%output(t)
END IF
END DO
END SUBROUTINE doProbes
END MODULE moduleProbe