diff --git a/doc/user-manual/.gitignore b/doc/user-manual/.gitignore index b2ba4d5..a8475d7 100644 --- a/doc/user-manual/.gitignore +++ b/doc/user-manual/.gitignore @@ -6,6 +6,7 @@ *.aux *.ps bibliography.bib.bak +bibliography.bib.sav *.bbl *.blg *.out diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index c424f50..1666f42 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index 2b8fa8e..e572331 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -72,7 +72,7 @@ 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. + The code is currently in the very early steps of development and further refinements are expected very soon. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Main Guidelines} @@ -86,11 +86,11 @@ \item \acrshort{fpakc} is coded in a \textit{understandable} way. This means that the code is required to be written in a clear way that is easy to understand and maintain. Variables and procedure names need to be self-understanding. - This ease the process of fixing bugs and improving the codes by a large team of developers. + This eases 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. + \item \acrshort{fpakc} requires being ease to use. 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. + \acrshort{fpakc} is aimed to be used in a wide range of applications and by various scientists: from well-established ones to newcomers to the field and also students. \end{enumerate} These are foundation stones of \acrshort{fpakc} and its development and should always be followed, at least for the releases in the official repository. @@ -105,16 +105,16 @@ \section{The Particle Method} \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. + For now own, macro-particles will be referred as just particles by abuse of language. 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. 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. + \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 are 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. @@ -124,10 +124,10 @@ \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. + \item If the number of iteration requires writing output files, it is done after all steps for the particles are completed. \end{enumerate} - \Gls{fpakc} has the capability to configure all the behavior of the simulation via the input file. + \Gls{fpakc} has the capability to configure all the behaviour 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}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -168,8 +168,8 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \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. + This is done by a neighbour search, starting from the previous cell containing the particle. + In the process of finding the new cell, a particle might encounter a boundary. 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. @@ -178,7 +178,7 @@ \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}. - These schemes detect when a particle change cells and split it if necessary to improve statistics. + These schemes detect when a particle changes cells and split it if necessary to improve statistics. The use of a Variable Weighting Scheme is defined by the user in the input file. Beware that this can increase the number of particles in the simulation and increase computational time. @@ -189,16 +189,16 @@ \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. + A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collision. \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 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 a random number is above the probability of collision for that specific type, the collision takes place. If not, the next type for the particle pair is checked. Below are described the type of collision process implemented in \acrshort{fpakc}: @@ -219,7 +219,7 @@ \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 modelled yet. + The photons emitted by this process are not modelled yet. \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -238,7 +238,7 @@ Once that the pushing is complete, the array of particles that remain inside the domain is copied to a new array. 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. + Unfortunately, this is done right now without parallelization and is very CPU consuming. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Probing}\label{sec:probing} @@ -250,21 +250,21 @@ 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. + If a particle velocity resides outside the velocity grid (in any direction), it will not be added to the tally of the distribution function. + Due to the limitation of only considering 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 considering 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} The properties of each particle are deposited in the nodes from the containing cell. - This process depend on the cell type, but in general, each node receive a proportional part of the particle properties as a function of the particle position inside the cell. - Figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes. + This process depends on the cell type, but in general, each node receives a proportional part of the particle properties as a function of the particle position inside the cell. + The figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes. \begin{wrapfigure}{l}{0.4\textwidth} \centering \includegraphics{figures/scatteringQuad} \caption{\label{fig:scatteringQuad}Example of how a particle is weighted in a quadrilateral cell.} \end{wrapfigure} - Each node receives a proportional part of the area formed by dividing the cell in for rectangles using as an additional vertex the particle position. + Each node receives a proportional part of the area formed by dividing the cell in for rectangles, using as an additional vertex the particle position. These properties are dimensionless, but they are converted to the correct units once the output is printed. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -273,11 +273,11 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Average scheme} - Particle-in-cell codes has an intrinsic statistical noise associated with them. + Particle-in-cell codes have an intrinsic statistical noise associated with them. Although this can be reduced by increasing the number of particles, this also increases the CPU requirements of the case. It is quite common that most cases reach a quasi-steady state after a number of iterations and time-average results can be obtained after to improve analysis, plotting and restarting the case using these time-average results as new species backgrounds. - Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the amount of iterations printed. + Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the number of iterations printed. \Gls{fpakc} implements a simple average scheme that, after a start time provided by the user, scores a mean and standard deviation of all the main species properties, and the electromagnetic field. This scheme is based on the Welford's online algorithm~\cite{welford1962note}. @@ -286,7 +286,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Installation} \section{Required Packages} - In order to properly compile \gls{fpakc}, the following packages are required. + To properly compile \gls{fpakc}, the following packages are required. \subsection{Gfortran} The \Gls{opensource} free compiler \Gls{gfortran}\cite{gfortranURL} from GCC is the basic way to compile \acrshort{fpakc}. It is distributed with all GNU/Linux distributions. @@ -369,7 +369,7 @@ make %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Case file} The required format for the case file is \Gls{json}. - \Gls{json} is a case-sensitive format, so input must be written with the correct capitalisation. + \Gls{json} is a case-sensitive format, so input must be written with the correct capitalization. The basic structure and options available for the case file are explained below. The order of the objects and variables is irrelevant, but the structure needs to be maintained. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -380,9 +380,9 @@ make \item \textbf{path}: Character. Path for the output files. This path is also used to locate the mesh input file. \item \textbf{folder}: Character. - Base name of the folder in wich output files are placed. + Base name of the folder in which output files are placed. The date and time is appended to this name. - If none is provided, only the date and time is writted as the folder name. + If none is provided, only the date and time is written as the folder name. \item \textbf{triggerOutput}: Integer. Determines the number of iterations between writing output files for macroscopic quantities. \item \textbf{cpuTime}: Logical. @@ -515,7 +515,7 @@ make \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. + \item \textbf{wallTemperature}: Reflective wall with constant temperature that exchange heat with particles. Required parameters are: \begin{itemize} \item \textbf{temperature}: Real. @@ -526,8 +526,8 @@ make 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. + A pair of ion-electron is generated for each ionization event, taking as a reference a neutral background. + The secondary electron is taken as the same type as the incident particle. The available input is: \begin{itemize} \item \textbf{neutral}: Object. @@ -540,7 +540,7 @@ make \item \textbf{mass}: Real. Units in $\unit{kg}$. Mass of neutral species. - If missing, the mass of the ion is ussed + If missing, the mass of the ion is used \item \textbf{density}: Real. Units in $\unit{m^{-3}}$. Density of neutral background. @@ -558,18 +558,18 @@ make \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. + As the particle is no longer simulated once it crossed the boundary, this time represents 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. + 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. + If , for some reason, a particle interacts with this axis, it is reflected. \end{itemize} \end{itemize} @@ -585,18 +585,26 @@ make Type of boundary. Accepted values are: \begin{itemize} - \item \textbf{dirichlet}: Elastic reflection of particles. + \item \textbf{dirichlet}: Constant value of electric potential on the surface. + \item \textbf{dirichletTime}: Constant value of the electric potential with a time variable profile. + The value of \textbf{boundaryEM.potential} will be multiplied for the corresponding value in the file \textbf{boundaryEM.temporalProfile}. \end{itemize} - \item \textbf{potential}: Real. - Fixed potential for Dirichlet boundary condition. + \item \textbf{potential}: Real. + Fixed potential for Dirichlet boundary condition. \item \textbf{physicalSurface}: Integer. Identification of the edge in the mesh file. + \item \textbf{temporalProfile}: Character. + Filename of the 2 column file containing the time variable profile. + File must be located in \textbf{output.path}. + The first column is the time in $\unit{s}$. + The second column is the factor that will multiply the value of the boundary. + \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{inject} The array \textbf{inject} specifies the injection of particles from different surfaces. - The injection of particles need to be associated to a physicalSurface in the mesh file. + The injection of particles needs to be associated to a physicalSurface in the mesh file. Multiple injections can be associated to the same surface. \begin{itemize} \item \textbf{name}: Character. @@ -610,7 +618,9 @@ make Available values are: \begin{itemize} \item \textbf{A}: Ampere. - \item \textbf{sccm}: Standard cubic centimeter. + \item \textbf{Am2}: Ampere per square meter. + This value will be multiplied by the area of injection. + \item \textbf{sccm}: Standard cubic centimetre. \item \textbf{part/s}: Particles (real) per second. \end{itemize} \item \textbf{v}: Real. @@ -627,7 +637,7 @@ make \begin{itemize} \item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction. \item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction. - Only takes into account the positive part of the half-Maxwellian. + Only considers the positive part of the half-Maxwellian. \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction. \end{itemize} \item \textbf{T}: Real. @@ -636,6 +646,11 @@ make Temperature in each direction. \item \textbf{physicalSurface}: Integer. Identification of the edge in the mesh file. + \item \textbf{particlesPerEdge}: Integer. + Optional. + Number of particles to be injected by each edge in the numerical domain. + The weight of the particles for each edge will modified by the surface of the edge to ensure the right flux is injected. + If no value is provided, the number of particles to inject per edge will be calculated with the species weight and the surface of the edge respect to the total one. \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{reference} @@ -651,7 +666,7 @@ make \item \textbf{radius}: Real. Reference atomic radius in $\unit{m}$. \item \textbf{crossSection}: Real. - Reference cross section in $\unit{m^2}$. + Reference cross-section in $\unit{m^2}$. If this value is present, radius is ignored. \end{itemize} @@ -677,8 +692,8 @@ make Indicates the type of pusher used for each species: \begin{itemize} \item \textbf{Neutral}: Pushes a particle without any external force. - \item \textbf{Electrostatic}: Pushes a particle including the effect of the electrostatic field. - \item \textbf{Electromagnetic}: Pushes particles accounting for the electromagnetic field. + \item \textbf{Electrostatic}: Pushes a particle, including the effect of the electrostatic field. + \item \textbf{Electromagnetic}: Pushes a particle, accounting for the electromagnetic field. \end{itemize} \item \textbf{WeightingScheme}: Character. Indicates the variable weighting scheme to be used in the simulation. @@ -706,11 +721,15 @@ make \begin{itemize} \item \textbf{species}: Character. Name of species as defined in the object \textbf{species}. - \item \textbf{file}: Character. - Output file from previous run used as an initial state for the species. - The file format must be the same as in \textbf{geometry.meshType} - Initial particles are assumed to have a Maxwellian distribution. - File must be located at \textbf{output.path}. + \item \textbf{file}: Character. + Output file from previous run used as an initial state for the species. + The file format must be the same as in \textbf{geometry.meshType} + Initial particles are assumed to have a Maxwellian distribution. + File must be located in \textbf{output.path}. + \item \textbf{particlesPerCell}: Integer. + Optional. + Initial number of particles per cell. + If not, the number of particles per cell will be assigned based on the species weight and the cell volume. \end{itemize} \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -726,11 +745,11 @@ make \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{interactions}\label{ssec:input_interactions} - This object determine the different interactions among species. + This object determines the different interactions among species. Acceptable values are: \begin{itemize} \item \textbf{folderCollisions}: Character. - Indicates the path to in which the cross section tables are allocated. + 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}. @@ -757,7 +776,7 @@ make Accepted values are \textbf{elastic}, \textbf{chargeExchange}, \textbf{ionization} and \textbf{recombination}. Please refer to Sec.~\ref{ssec:collisions} for a description of the different collision types. \item \textbf{crossSection}: Character. - File in \textbf{interactions.folderCollisions} that contains the cross section data as a 1D table of relative energy (in $\unit{eV}$) and cross section (in $\unit{m^-2}$). + File in \textbf{interactions.folderCollisions} that contains the cross-section data as a 1D table of relative energy (in $\unit{eV}$) and cross-section (in $\unit{m^-2}$). \item \textbf{energyThreshold}: Real. Energy threshold of the collisional process in $\unit{eV}$. Only valid for \textbf{ionization} and \textbf{recombination} processes. @@ -778,7 +797,7 @@ make \begin{itemize} \item \textbf{species\_i}, \textbf{species\_j}: Character. Define the two species involved in the collision processes. - Order is indiferent. + Order is indifferent. \end{itemize} \end{itemize} @@ -804,9 +823,9 @@ make %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{1D Emissive Cathode (1D\_Cathode)} - 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 modifying some parameters in the input file. + Emission from a 1D cathode in both, Cartesian and radial coordinates. + Both cases insert the same number of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field. + This case is useful to illustrate how \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. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/fpakc.f90 b/src/fpakc.f90 index e85f5f6..ae6e7bb 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -11,12 +11,12 @@ PROGRAM fpakc USE OMP_LIB IMPLICIT NONE - ! t = time step - INTEGER:: t ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 ! inputFile = path+name of input file CHARACTER(:), ALLOCATABLE:: inputFile + ! generic integer for time step + INTEGER:: t tStep = omp_get_wtime() !Gets the input file @@ -27,11 +27,18 @@ PROGRAM fpakc !Reads the json configuration file CALL readConfig(inputFile) + + !Create output folder and initial files + CALL initOutput(inputFile) + !Do '0' iteration - t = tInitial + timeStep = tInitial !$OMP PARALLEL DEFAULT(SHARED) !$OMP SINGLE + ! Initial reset of probes + CALL resetProbes() + CALL verboseError("Initial scatter of particles...") !$OMP END SINGLE CALL doScatter() @@ -45,19 +52,21 @@ PROGRAM fpakc tStep = omp_get_wtime() - tStep !Output initial state - CALL doOutput(t) + CALL doOutput() CALL verboseError('Starting main loop...') !$OMP PARALLEL DEFAULT(SHARED) DO t = tInitial + 1, tFinal - !Insert new particles and push them !$OMP SINGLE tStep = omp_get_wtime() + ! Update global time step index + timeStep = t + !Checks if a species needs to me moved in this iteration - CALL solver%updatePushSpecies(t) + CALL solver%updatePushSpecies() !Checks if probes need to be calculated this iteration - CALL resetProbes(t) + CALL resetProbes() tPush = omp_get_wtime() !$OMP END SINGLE @@ -75,7 +84,7 @@ PROGRAM fpakc !$OMP END SINGLE IF (doMCCollisions) THEN - CALL meshForMCC%doCollisions(t) + CALL meshForMCC%doCollisions() END IF @@ -120,12 +129,12 @@ PROGRAM fpakc !$OMP SINGLE tEMField = omp_get_wtime() - tEMField - CALL doAverage(t) + CALL doAverage() tStep = omp_get_wtime() - tStep !Output data - CALL doOutput(t) + CALL doOutput() !$OMP END SINGLE END DO diff --git a/src/modules/common/moduleCaseParam.f90 b/src/modules/common/moduleCaseParam.f90 index 8df3210..551d867 100644 --- a/src/modules/common/moduleCaseParam.f90 +++ b/src/modules/common/moduleCaseParam.f90 @@ -2,8 +2,13 @@ MODULE moduleCaseParam !Final and initial iterations INTEGER:: tFinal, tInitial = 0 + ! Global index of current iteration + INTEGER:: timeStep + ! Time step for all species REAL(8), ALLOCATABLE:: tau(:) + ! Minimum time step REAL(8):: tauMin + ! Time step for Monte-Carlo Collisions REAL(8):: tauColl END MODULE moduleCaseParam diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index cd553a8..ae5c548 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -40,10 +40,10 @@ MODULE moduleRandom INTEGER:: rnd REAL(8):: rnd01 - rnd = 0.D0 + rnd = 0 CALL RANDOM_NUMBER(rnd01) - rnd = INT(REAL(b - a) * rnd01) + 1 + rnd = a + FLOOR((b+1-a)*rnd01) END FUNCTION randomIntAB @@ -73,10 +73,21 @@ MODULE moduleRandom REAL(8), INTENT(in):: cumWeight(1:) REAL(8), INTENT(in):: sumWeight REAL(8):: rnd0b - INTEGER:: rnd + INTEGER:: rnd, i - rnd0b = random(0.D0, sumWeight) - rnd = MINLOC(DABS(rnd0b - cumWeight), 1) + rnd0b = random() + i = 1 + DO + IF (rnd0b <= cumWeight(i)/sumWeight) THEN + rnd = i + EXIT + + ELSE + i = i +1 + + END IF + END DO + ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1) END FUNCTION randomWeighted diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 7f6c036..5090942 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -84,20 +84,6 @@ MODULE moduleInput CALL readParallel(config) CALL checkStatus(config, "readParallel") - !If everything is correct, creates the output folder - CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder ) - !Copies input file to output folder - CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder) - !Copies particle mesh - IF (mesh%dimen > 0) THEN - CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder) - IF (doubleMesh) THEN - CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder) - - END IF - - END IF - END SUBROUTINE readConfig !Checks the status of the JSON case file and, if failed, exits the execution. @@ -282,8 +268,8 @@ MODULE moduleInput CALL readEMBoundary(config) !Read constant magnetic field DO i = 1, 3 - WRITE(istring, '(i2)') i - CALL config%get(object // '.B(' // istring // ')', B(i), found) + WRITE(iString, '(i2)') i + CALL config%get(object // '.B(' // iString // ')', B(i), found) IF (.NOT. found) THEN CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver') @@ -326,7 +312,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: object INTEGER:: nInitial - INTEGER:: i, j, p, e + INTEGER:: i, p, e CHARACTER(LEN=2):: iString CHARACTER(:), ALLOCATABLE:: spName INTEGER:: sp @@ -342,7 +328,8 @@ MODULE moduleInput REAL(8):: densityCen !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi - INTEGER:: nNewPart = 0.D0 + INTEGER:: nNewPart = 0 + REAL(8):: weight = 0.D0 CLASS(meshCell), POINTER:: cell TYPE(particle), POINTER:: partNew REAL(8):: vTh @@ -361,6 +348,9 @@ MODULE moduleInput !Reads node values at the nodes filename = path // spFile CALL mesh%readInitial(filename, density, velocity, temperature) + !Check if initial number of particles is given + CALL config%get(object // '.particlesPerCell', nNewPart, found) + !For each volume in the node, create corresponding particles DO e = 1, mesh%numCells !Scale variables @@ -373,7 +363,11 @@ MODULE moduleInput densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar) !Calculate number of particles - nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + IF (.NOT. found) THEN + nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + + END IF + weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart) !Allocate new particles DO p = 1, nNewPart @@ -410,7 +404,7 @@ MODULE moduleInput partNew%n_in = .TRUE. - partNew%weight = species(sp)%obj%weight + partNew%weight = weight !Assign particle to temporal list of particles CALL partInitial%add(partNew) @@ -809,7 +803,7 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config INTEGER:: i, s - CHARACTER(2):: istring, sString + CHARACTER(2):: iString, sString CHARACTER(:), ALLOCATABLE:: object, bType REAL(8):: Tw, cw !Wall temperature and specific heat !Neutral Properties @@ -817,16 +811,16 @@ MODULE moduleInput REAL(8), DIMENSION(:), ALLOCATABLE:: v0 REAL(8):: effTime REAL(8):: eThreshold !Energy threshold - INTEGER:: speciesID - CHARACTER(:), ALLOCATABLE:: speciesName, crossSection + INTEGER:: speciesID, electronSecondaryID + CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary LOGICAL:: found INTEGER:: nTypes CALL config%info('boundary', found, n_children = nBoundary) ALLOCATE(boundary(1:nBoundary)) DO i = 1, nBoundary - WRITE(istring, '(i2)') i - object = 'boundary(' // TRIM(istring) // ')' + WRITE(iString, '(i2)') i + object = 'boundary(' // TRIM(iString) // ')' boundary(i)%n = i CALL config%get(object // '.name', boundary(i)%name, found) @@ -873,8 +867,17 @@ MODULE moduleInput CALL config%get(object // '.crossSection', crossSection, found) IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') - CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & - speciesID, effTime, crossSection, eThreshold) + CALL config%get(object // '.electronSecondary', electronSecondary, found) + electronSecondaryID = speciesName2Index(electronSecondary) + IF (found) THEN + CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & + speciesID, effTime, crossSection, eThreshold,electronSecondaryID) + + ELSE + CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & + speciesID, effTime, crossSection, eThreshold) + + END IF CASE('wallTemperature') CALL config%get(object // '.temperature', Tw, found) @@ -924,7 +927,6 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile REAL(8):: volume - CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT object = 'geometry' @@ -1102,13 +1104,13 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(2):: istring + 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 + REAL(8):: everyTimeStep CALL config%info('output.probes', found, n_children = nProbes) @@ -1116,7 +1118,7 @@ MODULE moduleInput DO i = 1, nProbes WRITE(iString, '(I2)') i - object = 'output.probes(' // trim(istring) // ')' + object = 'output.probes(' // trim(iString) // ')' CALL config%get(object // '.species', speciesName, found) CALL config%get(object // '.position', r, found) @@ -1124,16 +1126,14 @@ MODULE moduleInput 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) + CALL config%get(object // '.timeStep', everyTimeStep, 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) + CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep) END DO - CALL resetProbes(tInitial) - END SUBROUTINE readProbes SUBROUTINE readEMBoundary(config) @@ -1141,7 +1141,6 @@ MODULE moduleInput USE moduleOutput USE moduleErrors USE moduleEM - USE moduleRefParam USE moduleSpecies USE json_module IMPLICIT NONE @@ -1149,34 +1148,72 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(2):: istring - INTEGER:: i, e, s + CHARACTER(:), ALLOCATABLE:: typeEM + REAL(8):: potential + INTEGER:: physicalSurface + CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath + INTEGER:: b, s, n, ni + CHARACTER(2):: bString INTEGER:: info EXTERNAL:: dgetrf CALL config%info('boundaryEM', found, n_children = nBoundaryEM) - IF (found) ALLOCATE(boundEM(1:nBoundaryEM)) + IF (found) THEN + ALLOCATE(boundaryEM(1:nBoundaryEM)) + + END IF - DO i = 1, nBoundaryEM - WRITE(istring, '(I2)') i - object = 'boundaryEM(' // trim(istring) // ')' + DO b = 1, nBoundaryEM + WRITE(bString, '(I2)') b + object = 'boundaryEM(' // TRIM(bString) // ')' - CALL config%get(object // '.type', boundEM(i)%typeEM, found) + CALL config%get(object // '.type', typeEM, found) - SELECT CASE(boundEM(i)%typeEM) + SELECT CASE(typeEM) CASE ("dirichlet") - CALL config%get(object // '.potential', boundEM(i)%potential, found) - IF (.NOT. found) & + CALL config%get(object // '.potential', potential, found) + IF (.NOT. found) THEN CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') - boundEM(i)%potential = boundEM(i)%potential/Volt_ref - CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found) - IF (.NOT. found) & - CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') + END IF + + CALL config%get(object // '.physicalSurface', physicalSurface, found) + IF (.NOT. found) THEN + CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', & + 'readEMBoundary') + + END IF + + CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential) + + CASE ("dirichletTime") + CALL config%get(object // '.potential', potential, found) + IF (.NOT. found) THEN + CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', & + 'readEMBoundary') + + END IF + + CALL config%get(object // '.temporalProfile', temporalProfile, found) + IF (.NOT. found) THEN + CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', & + 'readEMBoundary') + + END IF + temporalProfilePath = path // temporalProfile + + CALL config%get(object // '.physicalSurface', physicalSurface, found) + IF (.NOT. found) THEN + CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', & + 'readEMBoundary') + + END IF + + CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath) CASE DEFAULT - CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary') + CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') END SELECT @@ -1195,18 +1232,28 @@ MODULE moduleInput END DO - IF (ALLOCATED(boundEM)) THEN - DO e = 1, mesh%numEdges - IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN - DO i = 1, nBoundaryEM - IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN - CALL boundEM(i)%apply(mesh%edges(e)%obj) + ! Modify K matrix due to boundary conditions + DO b = 1, nBoundaryEM + SELECT TYPE(boundary => boundaryEM(b)%obj) + TYPE IS(boundaryEMDirichlet) + DO n = 1, boundary%nNodes + ni = boundary%nodes(n)%obj%n + mesh%K(ni, :) = 0.D0 + mesh%K(ni, ni) = 1.D0 - END IF - END DO - END IF - END DO - END IF + END DO + + TYPE IS(boundaryEMDirichletTime) + DO n = 1, boundary%nNodes + ni = boundary%nodes(n)%obj%n + mesh%K(ni, :) = 0.D0 + mesh%K(ni, ni) = 1.D0 + + END DO + + END SELECT + + END DO !Compute the PLU factorization of K once boundary conditions have been read CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info) @@ -1227,24 +1274,25 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config INTEGER:: i - CHARACTER(2):: istring + CHARACTER(2):: iString CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found CHARACTER(:), ALLOCATABLE:: speciesName CHARACTER(:), ALLOCATABLE:: name REAL(8):: v - REAL(8), ALLOCATABLE:: T(:), normal(:) + REAL(8), ALLOCATABLE:: temperature(:), normal(:) REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface + INTEGER:: particlesPerEdge INTEGER:: sp CALL config%info('inject', found, n_children = nInject) ALLOCATE(inject(1:nInject)) nPartInj = 0 DO i = 1, nInject - WRITE(istring, '(i2)') i - object = 'inject(' // trim(istring) // ')' + WRITE(iString, '(i2)') i + object = 'inject(' // trim(iString) // ')' !Find species CALL config%get(object // '.species', speciesName, found) @@ -1252,7 +1300,7 @@ MODULE moduleInput CALL config%get(object // '.name', name, found) CALL config%get(object // '.v', v, found) - CALL config%get(object // '.T', T, found) + CALL config%get(object // '.T', temperature, found) CALL config%get(object // '.n', normal, found) IF (.NOT. found) THEN ALLOCATE(normal(1:3)) @@ -1261,8 +1309,10 @@ MODULE moduleInput CALL config%get(object // '.flow', flow, found) CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) + particlesPerEdge = 0 + CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) - CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) + CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge) CALL readVelDistr(config, inject(i), object) @@ -1322,28 +1372,28 @@ MODULE moduleInput TYPE(injectGeneric), INTENT(inout):: inj CHARACTER(:), ALLOCATABLE, INTENT(in):: object INTEGER:: i - CHARACTER(2):: istring + CHARACTER(2):: iString CHARACTER(:), ALLOCATABLE:: fvType LOGICAL:: found - REAL(8):: v, T, m + REAL(8):: v, temperature, m !Reads species mass m = inj%species%m !Reads distribution functions for velocity DO i = 1, 3 - WRITE(istring, '(i2)') i - CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found) - IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // & + WRITE(iString, '(i2)') i + CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found) + IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // & " found for " // object, 'readVelDistr') SELECT CASE(fvType) CASE ("Maxwellian") - T = inj%T(i) - CALL initVelDistMaxwellian(inj%v(i)%obj, t, m) + temperature = inj%temperature(i) + CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m) CASE ("Half-Maxwellian") - T = inj%T(i) - CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m) + temperature = inj%temperature(i) + CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m) CASE ("Delta") v = inj%vMod*inj%n(i) @@ -1384,5 +1434,37 @@ MODULE moduleInput END SUBROUTINE readParallel + SUBROUTINE initOutput(inputFile) + USE moduleRefParam + USE moduleMesh, ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl + USE moduleOutput, ONLY: path, folder + IMPLICIT NONE + + CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile + INTEGER:: fileReference = 30 + !If everything is correct, creates the output folder + CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder ) + !Copies input file to output folder + CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder) + !Copies particle mesh + IF (mesh%dimen > 0) THEN + CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder) + IF (doubleMesh) THEN + CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder) + + END IF + + END IF + + ! Write commit of fpakc + CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpakc_commit.txt') + + ! Write file with reference values + OPEN (fileReference, file=path // folder // '/' // 'reference.txt') + WRITE(fileReference, "(7(1X,A20))") 'L_ref', 'v_ref', 'ti_ref', 'Vol_ref', 'EF_ref', 'Volt_ref', 'B_ref' + WRITE(fileReference, "(7(1X,ES20.6E3))") L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref + CLOSE(fileReference) + + END SUBROUTINE initOutput END MODULE moduleInput diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 269f157..bbc72e2 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -104,6 +104,7 @@ MODULE moduleMesh1DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(out):: self @@ -122,6 +123,8 @@ MODULE moduleMesh1DCart self%x = r1(1) + self%surface = 1.D0 + self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index d998267..e260900 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -104,6 +104,7 @@ MODULE moduleMesh1DRad USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(out):: self @@ -122,6 +123,8 @@ MODULE moduleMesh1DRad self%r = r1(1) + self%surface = 1.D0 + self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 4341cb0..dbc8b25 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -144,6 +144,7 @@ MODULE moduleMesh2DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(out):: self @@ -163,7 +164,7 @@ MODULE moduleMesh2DCart r2 = self%n2%getCoordinates() self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) - self%weight = 1.D0 + self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -318,6 +319,8 @@ MODULE moduleMesh2DCart INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) + fPsi = 0.D0 + fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & @@ -508,15 +511,15 @@ MODULE moduleMesh2DCart conv = 1.D0 XiO = 0.D0 + f(3) = 0.D0 DO WHILE(conv > 1.D-4) dPsi = self%dPsi(XiO, 4) pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) fPsi = self%fPsi(XiO, 4) - f = (/ DOT_PRODUCT(fPsi,self%x), & - DOT_PRODUCT(fPsi,self%y), & - 0.D0 /) - r + f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & + DOT_PRODUCT(fPsi,self%y) /) - r(1:2) Xi = XiO - MATMUL(invJ, f)/detJ conv = MAXVAL(DABS(Xi-XiO),1) XiO = Xi @@ -554,6 +557,7 @@ MODULE moduleMesh2DCart !Compute element volume PURE SUBROUTINE volumeQuad(self) + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(inout):: self @@ -569,8 +573,9 @@ MODULE moduleMesh2DCart pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) + !Compute total volume of the cell - self%volume = detJ*4.D0 + self%volume = detJ*4.D0/L_ref !Compute volume per node self%n1%v = self%n1%v + fPsi(1)*self%volume self%n2%v = self%n2%v + fPsi(2)*self%volume @@ -762,6 +767,7 @@ MODULE moduleMesh2DCart pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index d4baedd..ae1eb92 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -144,6 +144,7 @@ MODULE moduleMesh2DCyl USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleConstParam, ONLY: PI IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(out):: self @@ -163,7 +164,15 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = r2(2)**2 - r1(2)**2 + !Edge surface + IF (self%z(2) /= self%z(1)) THEN + self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1)) + + ELSE + self%surface = ABS(self%r(2)**2 - self%r(1)**2) + + END IF + self%surface = self%surface * PI !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & @@ -223,21 +232,13 @@ MODULE moduleMesh2DCyl CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd REAL(8):: r(1:3) - REAL(8):: dr, dz + REAL(8):: p1(1:2), p2(1:2) rnd = random() - dr = self%r(2) - self%r(1) - dz = self%z(2) - self%z(1) - IF (dr /= 0.D0) THEN - r(2) = dr * DSQRT(rnd) + self%r(1) - r(1) = dz * (r(2) - self%r(1))/dr + self%z(1) - - ELSE - r(2) = self%r(1) - r(1) = dz * rnd + self%z(1) - - END IF + p1 = (/self%z(1), self%r(1) /) + p2 = (/self%z(2), self%r(2) /) + r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 r(3) = 0.D0 END FUNCTION randPosEdge @@ -246,7 +247,6 @@ MODULE moduleMesh2DCyl !QUAD FUNCTIONS !Init element SUBROUTINE initCellQuad2DCyl(self, n, p, nodes) - USE moduleRefParam IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(out):: self @@ -326,6 +326,8 @@ MODULE moduleMesh2DCyl INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) + fPsi = 0.D0 + fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & @@ -496,7 +498,7 @@ MODULE moduleMesh2DCyl END FUNCTION elemFQuad - !Checks if Xi is inside the element + !Check if Xi is inside the element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -524,15 +526,15 @@ MODULE moduleMesh2DCyl conv = 1.D0 XiO = 0.D0 + f(3) = 0.D0 DO WHILE(conv > 1.D-4) dPsi = self%dPsi(XiO, 4) pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) fPsi = self%fPsi(XiO, 4) - f = (/ DOT_PRODUCT(fPsi,self%z), & - DOT_PRODUCT(fPsi,self%r), & - 0.D0 /) - r + f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), & + DOT_PRODUCT(fPsi,self%r) /) - r(1:2) Xi = XiO - MATMUL(invJ, f)/detJ conv = MAXVAL(DABS(Xi-XiO),1) XiO = Xi @@ -553,7 +555,7 @@ MODULE moduleMesh2DCyl XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) nextInt = MAXLOC(XiArray,1) - !Selects the higher value of directions and searches in that direction + !Select the higher value of directions and searches in that direction NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) @@ -581,6 +583,7 @@ MODULE moduleMesh2DCyl REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) self%volume = 0.D0 + !2D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 4) @@ -589,18 +592,18 @@ MODULE moduleMesh2DCyl fPsi = self%fPsi(Xi, 4) r = DOT_PRODUCT(fPsi,self%r) !Computes total volume of the cell - self%volume = r*detJ*PI8 !4*2*pi - !Computes volume per node - Xi = (/-5.D-1, -5.D-1, 0.D0/) + self%volume = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral) + !Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis. + Xi = (/-5.D-1, -0.25D0, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 - Xi = (/ 5.D-1, -5.D-1, 0.D0/) + Xi = (/ 5.D-1, -0.25D0, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 - Xi = (/ 5.D-1, 5.D-1, 0.D0/) + Xi = (/ 5.D-1, 0.75D0, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 - Xi = (/-5.D-1, 5.D-1, 0.D0/) + Xi = (/-5.D-1, 0.75D0, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index c451689..8170df7 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -109,6 +109,7 @@ MODULE moduleMesh3DCart USE moduleBoundary USE moduleErrors USE moduleMath + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(out):: self @@ -142,6 +143,8 @@ MODULE moduleMesh3DCart self%normal = crossProduct(vec1, vec2) self%normal = normalize(self%normal) + self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D + !Boundary index self%boundary => boundary(bt) ALLOCATE(self%fBoundary(1:nSpecies)) diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index 97ec729..c0dcfbb 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -1,22 +1,22 @@ MODULE moduleMeshOutput0D CONTAINS - SUBROUTINE printOutput0D(self, t) + SUBROUTINE printOutput0D(self) USE moduleMesh USE moduleRefParam USE moduleSpecies USE moduleOutput + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: i TYPE(outputFormat):: output CHARACTER(:), ALLOCATABLE:: fileName DO i = 1, nSpecies fileName='OUTPUT_' // species(i)%obj%name // '.dat' - IF (t == 0) THEN + IF (timeStep == 0) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", & "pressure (Pa)", "temperature (K)" @@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj) - WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature + WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density, & + output%velocity, & + output%pressure, & + output%temperature CLOSE(20) END DO END SUBROUTINE printOutput0D - SUBROUTINE printColl0D(self, t) + SUBROUTINE printColl0D(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName INTEGER:: k fileName='OUTPUT_Collisions.dat' - IF (t == tInitial) THEN + IF (timeStep == tInitial) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions" WRITE(*, "(6X,A15,A)") "Creating file: ", fileName @@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) + WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) CLOSE(20) END SUBROUTINE printColl0D - SUBROUTINE printEM0D(self, t) + SUBROUTINE printEM0D(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printEM0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 0df1289..1ad6a36 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -108,6 +108,7 @@ MODULE moduleMeshInputGmsh2 READ(10, *) totalNumElem !Count edges and volume elements + numEdges = 0 SELECT TYPE(self) TYPE IS(meshParticles) self%numEdges = 0 @@ -328,7 +329,7 @@ MODULE moduleMeshInputGmsh2 DO i = 1, numNodes !Reads the density - READ(10, *), e, density(i) + READ(10, *) e, density(i) END DO @@ -339,7 +340,7 @@ MODULE moduleMeshInputGmsh2 DO i = 1, numNodes !Reads the velocity - READ(10, *), e, velocity(i, 1:3) + READ(10, *) e, velocity(i, 1:3) END DO diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index ccdcf3d..8176bb5 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE writeGmsh2FooterElementData !Prints the scattered properties of particles into the nodes - SUBROUTINE printOutputGmsh2(self, t) + SUBROUTINE printOutputGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleSpecies USE moduleOutput USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: n, i TYPE(outputFormat):: output(1:self%numNodes) REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t) + fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) WRITE(60, "(I6,ES20.6E3)") n, output(n)%density END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature END DO @@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printOutputGmsh2 !Prints the number of collisions into the volumes - SUBROUTINE printCollGmsh2(self, t) + SUBROUTINE printCollGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: k, c INTEGER:: n @@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2 END SELECT IF (collOutput) THEN - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref - fileName = formatFileName(prefix, 'Collisions', 'msh', t) + fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2 DO c = 1, interactionMatrix(k)%amount WRITE(cString, "(I2)") c title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString - CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells) + CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells) DO n=1, self%numCells WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c) END DO @@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printCollGmsh2 !Prints the electrostatic EM properties into the nodes and volumes - SUBROUTINE printEMGmsh2(self, t) + SUBROUTINE printEMGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: n, e REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName @@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2 Xi = (/ 0.D0, 0.D0, 0.D0 /) IF (emOutput) THEN - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref - fileName = formatFileName(prefix, 'EMField', 'msh', t) + fileName = formatFileName(prefix, 'EMField', 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (20, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(20) - CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref END DO CALL writeGmsh2FooterNodeData(20) - CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells) + CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells) DO e=1, self%numCells WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref END DO CALL writeGmsh2FooterElementData(20) - CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref END DO diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 index e4a6c72..7dc2e84 100644 --- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 +++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 @@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon CHARACTER(LEN=4):: prefix = 'Step' CONTAINS - PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName) + PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName) USE moduleOutput IMPLICIT NONE CHARACTER(*), INTENT(in):: prefix, suffix, extension - INTEGER, INTENT(in), OPTIONAL:: t + INTEGER, INTENT(in), OPTIONAL:: timeStep CHARACTER (LEN=iterationDigits):: tString CHARACTER(:), ALLOCATABLE:: fileName - IF (PRESENT(t)) THEN - WRITE(tString, iterationFormat) t + IF (PRESENT(timeStep)) THEN + WRITE(tString, iterationFormat) timeStep fileName = prefix // '_' // tString // '_' // suffix // '.' // extension ELSE diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 763517f..e07db01 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU CLASS(meshGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8):: r(1:3) !3 generic coordinates - INTEGER:: fileID, error, found + INTEGER:: fileID CHARACTER(LEN=256):: line INTEGER:: numNodes, numElements, numEdges INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types @@ -275,6 +275,7 @@ MODULE moduleMeshInputVTU END DO !Count the number of edges + numEdges = 0 SELECT CASE(self%dimen) CASE(3) !Edges are triangles, type 5 in VTK @@ -495,7 +496,7 @@ MODULE moduleMeshInputVTU END SELECT END DO - + !Call mesh connectivity CALL self%connectMesh diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 index 6286cfc..81e4bbf 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -209,23 +209,22 @@ MODULE moduleMeshOutputVTU WRITE(fileID,"(8X,A)") '' !Electric field WRITE(fileID,"(10X,A, A, A)") '' - WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' END SUBROUTINE writeEM - SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection) + SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection) USE moduleCaseParam USE moduleOutput USE moduleRefParam IMPLICIT NONE INTEGER:: fileID - INTEGER, INTENT(in):: t CHARACTER(*):: fileNameStep, fileNameCollection - IF (t == tInitial) THEN + IF (timeStep == tInitial) THEN !Create collection file WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection) @@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU !Write iteration file in collection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND') - WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '' + WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") & + '' !Close collection file - IF (t == tFinal) THEN + IF (timeStep == tFinal) THEN WRITE (fileID + 1, "(2X, A)") '' WRITE (fileID + 1, "(A)") '' @@ -307,22 +307,21 @@ MODULE moduleMeshOutputVTU END SUBROUTINE writeAverage - SUBROUTINE printOutputVTU(self,t) + SUBROUTINE printOutputVTU(self) USE moduleMesh USE moduleSpecies USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n, i, fileID + INTEGER:: i, fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection - TYPE(outputFormat):: output(1:self%numNodes) fileID = 60 DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t) + fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -338,29 +337,27 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END DO END SUBROUTINE printOutputVTU - SUBROUTINE printCollVTU(self,t) + SUBROUTINE printCollVTU(self) USE moduleMesh USE moduleOutput USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n, i, fileID + INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection - CHARACTER (LEN=iterationDigits):: tstring - TYPE(outputFormat):: output(1:self%numNodes) fileID = 62 IF (collOutput) THEN - fileName = formatFileName(prefix, 'Collisions', 'vtu', t) + fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -376,26 +373,26 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END IF END SUBROUTINE printCollVTU - SUBROUTINE printEMVTU(self, t) + SUBROUTINE printEMVTU(self) USE moduleMesh USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection fileID = 64 IF (emOutput) THEN - fileName = formatFileName(prefix, 'EMField', 'vtu', t) + fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -411,7 +408,7 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'EMField', 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END IF @@ -424,9 +421,8 @@ MODULE moduleMeshOutputVTU IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER:: n, i, fileIDMean, fileIDDeviation + INTEGER:: i, fileIDMean, fileIDDeviation CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation - TYPE(outputFormat):: output(1:self%numNodes) fileIDMean = 66 fileIDDeviation = 67 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index e96ff2a..7ab3914 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -59,6 +59,13 @@ MODULE moduleMesh END TYPE meshNodeCont + ! Array of pointers to nodes. + TYPE:: meshNodePointer + CLASS(meshNode), POINTER:: obj + CONTAINS + + END TYPE meshNodePointer + !Type for array of boundary functions (one per species) TYPE, PUBLIC:: fBoundaryGeneric PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL() @@ -76,8 +83,8 @@ MODULE moduleMesh CLASS(meshCell), POINTER:: eColl => NULL() !Normal vector REAL(8):: normal(1:3) - !Weight for random injection of particles - REAL(8):: weight = 1.D0 + ! Surface of edge + REAL(8):: surface = 0.D0 !Pointer to boundary type TYPE(boundaryCont), POINTER:: boundary !Array of functions for boundary conditions @@ -372,10 +379,9 @@ MODULE moduleMesh END SUBROUTINE connectMesh_interface !Prints number of collisions in each cell - SUBROUTINE printColl_interface(self, t) + SUBROUTINE printColl_interface(self) IMPORT meshGeneric CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface @@ -403,18 +409,16 @@ MODULE moduleMesh ABSTRACT INTERFACE !Prints Species data - SUBROUTINE printOutput_interface(self, t) + SUBROUTINE printOutput_interface(self) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printOutput_interface !Prints EM info - SUBROUTINE printEM_interface(self, t) + SUBROUTINE printEM_interface(self) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface @@ -613,6 +617,7 @@ MODULE moduleMesh INTEGER:: sp INTEGER:: i CLASS(meshNode), POINTER:: node + REAL(8):: pFraction !Particle fraction cellNodes = self%getNodes(nNodes) fPsi = self%fPsi(part%Xi, nNodes) @@ -623,10 +628,11 @@ MODULE moduleMesh DO i = 1, nNodes node => mesh%nodes(cellNodes(i))%obj + pFraction = fPsi(i)*part%weight CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS + node%output(sp)%den = node%output(sp)%den + pFraction + node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:) + node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS CALL OMP_UNSET_LOCK(node%lock) END DO @@ -787,7 +793,7 @@ MODULE moduleMesh END FUNCTION findCellBrute !Computes collisions in element - SUBROUTINE doCollisions(self, t) + SUBROUTINE doCollisions(self) USE moduleCollisions USE moduleSpecies USE moduleList @@ -795,10 +801,10 @@ MODULE moduleMesh USE moduleRandom USE moduleOutput USE moduleMath + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(inout), TARGET:: self - INTEGER, INTENT(in):: t INTEGER:: e CLASS(meshCell), POINTER:: cell INTEGER:: k, i, j @@ -814,7 +820,7 @@ MODULE moduleMesh REAL(8):: rnd_real !Random number for collision INTEGER:: rnd_int !Random number for collision - IF (MOD(t, everyColl) == 0) THEN + IF (MOD(timeStep, everyColl) == 0) THEN !Collisions need to be performed in this iteration !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) DO e=1, self%numCells @@ -1023,6 +1029,9 @@ MODULE moduleMesh ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) + deltaV_ij = 0.D0 + p_ij = 0.D0 + mass_ij = 0.D0 !Loop over particles of species_i partTemp => cell%listPart_in(i)%head p = 1 @@ -1107,6 +1116,9 @@ MODULE moduleMesh ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) + deltaV_ji = 0.D0 + p_ji = 0.D0 + mass_ji = 0.D0 !Loop over particles of species_j partTemp => cell%listPart_in(j)%head p = 1 diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 4f72e10..091e52e 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -147,7 +147,13 @@ MODULE moduleMeshBoundary ALLOCATE(newElectron) ALLOCATE(newIon) - newElectron%species => part%species + IF (ASSOCIATED(bound%electronSecondary)) THEN + newElectron%species => bound%electronSecondary + + ELSE + newElectron%species => part%species + + END IF newIon%species => bound%species newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0)) diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 83c815c..0b76105 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -38,6 +38,7 @@ MODULE moduleBoundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. CLASS(speciesGeneric), POINTER:: species !Ion species + CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron TYPE(table1D):: crossSection REAL(8):: effectiveTime REAL(8):: eThreshold @@ -103,17 +104,19 @@ MODULE moduleBoundary END SUBROUTINE initWallTemperature - SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold) + SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) USE moduleRefParam USE moduleSpecies USE moduleCaseParam USE moduleConstParam + USE moduleErrors IMPLICIT NONE CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary REAL(8), INTENT(in):: me !Electron mass REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties - INTEGER:: speciesID + INTEGER, INTENT(in):: ion + INTEGER, OPTIONAL, INTENT(in):: electronSecondary REAL(8):: effTime CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection REAL(8), INTENT(in):: eThreshold @@ -126,7 +129,22 @@ MODULE moduleBoundary boundary%n0 = n0 * Vol_ref boundary%v0 = v0 / v_ref boundary%vTh = DSQRT(kb*T0/m0)/v_ref - boundary%species => species(speciesID)%obj + boundary%species => species(ion)%obj + IF (PRESENT(electronSecondary)) THEN + SELECT TYPE(sp => species(electronSecondary)%obj) + TYPE IS(speciesCharged) + boundary%electronSecondary => sp + + CLASS DEFAULT + CALL criticalError("Species " // sp%name // " chosen for " // & + "secondary electron is not a charged species", 'initIonization') + + END SELECT + + ELSE + boundary%electronSecondary => NULL() + + END IF boundary%effectiveTime = effTime / ti_ref CALL boundary%crossSection%init(crossSection) CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 496ea6a..b4e0ed2 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -54,15 +54,16 @@ MODULE moduleInject INTEGER:: id CHARACTER(:), ALLOCATABLE:: name REAL(8):: vMod !Velocity (module) - REAL(8):: T(1:3) !Temperature + REAL(8):: temperature(1:3) !Temperature REAL(8):: n(1:3) !Direction of injection LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n INTEGER:: nParticles !Number of particles to introduce each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges - REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability - REAL(8):: sumWeight + INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge + REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge + REAL(8):: surface ! Total surface of injection TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction CONTAINS PROCEDURE, PASS:: init => initInject @@ -75,7 +76,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) + SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -86,50 +87,29 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i - REAL(8), INTENT(in):: v, n(1:3), T(1:3) - INTEGER, INTENT(in):: sp, physicalSurface + REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) + INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge REAL(8):: tauInject REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) INTEGER:: nVolColl + REAL(8):: fluxPerStep = 0.D0 - self%id = i - self%vMod = v / v_ref - self%n = n / NORM2(n) - self%T = T / T_ref - self%species => species(sp)%obj - tauInject = tau(self%species%n) - SELECT CASE(units) - CASE ("sccm") - !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight) - - CASE ("A") - !Input current in Ampers - ! TODO: Make this only available for charge species - self%nParticles = INT(flow*tauInject*ti_ref/(qe*abs(species(sp)%obj%qm*species(sp)%obj%m)*species(sp)%obj%weight)) - - CASE ("part/s") - !Input current in Ampers - self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight) - - CASE DEFAULT - CALL criticalError("No support for units: " // units, 'initInject') - - END SELECT - !Scale particles for different species steps - IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') - + self%id = i + self%vMod = v / v_ref + self%n = n / NORM2(n) + self%temperature = temperature / T_ref !Gets the edge elements from which particles are injected DO e = 1, mesh%numEdges phSurface(e) = mesh%edges(e)%obj%physicalSurface END DO - self%nEdges = COUNT(phSurface == physicalSurface) - ALLOCATE(inject(i)%edges(1:self%nEdges)) + ALLOCATE(self%edges(1:self%nEdges)) + ALLOCATE(self%particlesPerEdge(1:self%nEdges)) + ALLOCATE(self%weightPerEdge(1:self%nEdges)) et = 0 DO e=1, mesh%numEdges IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN @@ -161,15 +141,78 @@ MODULE moduleInject END DO - !Calculates cumulative probability - ALLOCATE(self%cumWeight(1:self%nEdges)) - et = 1 - self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight - DO et = 2, self%nEdges - self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1) + !Calculates total area + self%surface = 0.D0 + DO et = 1, self%nEdges + self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface END DO - self%sumWeight = self%cumWeight(self%nEdges) + + ! Information about species and flux + self%species => species(sp)%obj + tauInject = tau(self%species%n) + ! Convert units + SELECT CASE(units) + CASE ("sccm") + !Standard cubic centimeter per minute + fluxPerStep = flow*sccm2atomPerS + + CASE ("A") + !Current in Ampers + SELECT TYPE(sp => self%species) + CLASS IS(speciesCharged) + fluxPerStep = flow/(qe*abs(sp%q)) + + CLASS DEFAULT + call criticalError('Attempted to assign a flux in "A" to a species without charge.', 'initInject') + + END SELECT + + CASE ("Am2") + !Input current in Ampers per square meter + SELECT TYPE(sp => self%species) + CLASS IS(speciesCharged) + fluxPerStep = flow*self%surface*L_ref**2/(qe*abs(sp%q)) + + CLASS DEFAULT + call criticalError('Attempted to assign a flux in "Am2" to a species without charge.', 'initInject') + + END SELECT + + + CASE ("part/s") + !Input current in Ampers + fluxPerStep = flow + + CASE DEFAULT + CALL criticalError("No support for units: " // units, 'initInject') + + END SELECT + fluxPerStep = fluxPerStep * tauInject * ti_ref / self%surface + + !Assign particles per edge + IF (particlesPerEdge > 0) THEN + ! Particles per edge defined by the user + self%particlesPerEdge = particlesPerEdge + DO et = 1, self%nEdges + self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge) + + END DO + + ELSE + ! No particles assigned per edge, use the species weight + self%weightPerEdge = self%species%weight + DO et = 1, self%nEdges + self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight) + + END DO + + END IF + + self%nParticles = SUM(self%particlesPerEdge) + + !Scale particles for different species steps + IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') END SUBROUTINE initInject @@ -204,23 +247,23 @@ MODULE moduleInject END SUBROUTINE doInjects - SUBROUTINE initVelDistMaxwellian(velDist, T, m) + SUBROUTINE initVelDistMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: T, m + REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(T/m)) + velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistMaxwellian - SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m) + SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: T, m + REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -283,9 +326,8 @@ MODULE moduleInject IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - INTEGER:: randomX - INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array - INTEGER:: i + INTEGER, SAVE:: nMin + INTEGER:: i, e INTEGER:: n, sp CLASS(meshEdge), POINTER:: randomEdge REAL(8):: direction(1:3) @@ -300,59 +342,62 @@ MODULE moduleInject END IF END DO - nMin = nMin + 1 - nMax = nMin + self%nParticles - 1 - !Assign weight to particle. - partInj(nMin:nMax)%weight = self%species%weight - !Particle is considered to be outside the domain - partInj(nMin:nMax)%n_in = .FALSE. !$OMP END SINGLE !$OMP DO - DO n = nMin, nMax - randomX = randomWeighted(self%cumWeight, self%sumWeight) + DO e = 1, self%nEdges + ! Select edge for injection + randomEdge => mesh%edges(self%edges(e))%obj + ! Inject particles in edge + DO i = 1, self%particlesPerEdge(e) + ! Index in the global partInj array + n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i + !Particle is considered to be outside the domain + partInj(n)%n_in = .FALSE. + !Random position in edge + partInj(n)%r = randomEdge%randPos() + !Assign weight to particle. + partInj(n)%weight = self%weightPerEdge(e) + !Volume associated to the edge: + IF (ASSOCIATED(randomEdge%e1)) THEN + partInj(n)%cell = randomEdge%e1%n - randomEdge => mesh%edges(self%edges(randomX))%obj - !Random position in edge - partInj(n)%r = randomEdge%randPos() - !Volume associated to the edge: - IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%cell = randomEdge%e1%n + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN + partInj(n)%cell = randomEdge%e2%n - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n + ELSE + CALL criticalError("No Volume associated to edge", 'addParticles') - ELSE - CALL criticalError("No Volume associated to edge", 'addParticles') + END IF + partInj(n)%cellColl = randomEdge%eColl%n + sp = self%species%n - END IF - partInj(n)%cellColl = randomEdge%eColl%n - sp = self%species%n + !Assign particle type + partInj(n)%species => self%species - !Assign particle type - partInj(n)%species => self%species + direction = self%n - direction = self%n + partInj(n)%v = 0.D0 - partInj(n)%v = 0.D0 + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) + !If velocity is not in the right direction, invert it + IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN + partInj(n)%v = - partInj(n)%v - !If velocity is not in the right direction, invert it - IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN - partInj(n)%v = - partInj(n)%v + END IF - END IF + !Obtain natural coordinates of particle in cell + partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) + !Push new particle with the minimum time step + CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) + !Assign cell to new particle + CALL solver%updateParticleCell(partInj(n)) - !Obtain natural coordinates of particle in cell - partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) - !Push new particle with the minimum time step - CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) - !Assign cell to new particle - CALL solver%updateParticleCell(partInj(n)) + END DO END DO !$OMP END DO diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90 index c7d3cf5..754d56a 100644 --- a/src/modules/moduleProbe.f90 +++ b/src/modules/moduleProbe.f90 @@ -27,7 +27,7 @@ MODULE moduleProbe CONTAINS !Functions for probeDistFunc type - SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep) + SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep) USE moduleCaseParam USE moduleRefParam USE moduleSpecies @@ -41,7 +41,7 @@ MODULE moduleProbe 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 + REAL(8), INTENT(in):: everyTimeStep INTEGER:: sp, i REAL(8):: dv(1:3) @@ -91,17 +91,17 @@ MODULE moduleProbe 1:self%nv(3))) !Number of iterations between output - IF (timeStep == 0.D0) THEN + IF (everyTimeStep == 0.D0) THEN self%every = 1 ELSE - self%every = NINT(timeStep/ tauMin / ti_ref) + self%every = NINT(everyTimeStep/ tauMin / ti_ref) END IF !Maximum radius !TODO: Make this an input parameter - self%maxR = 1.D0 + self%maxR = 1.D-2/L_ref !Init the probe lock CALL OMP_INIT_LOCK(self%lock) @@ -148,7 +148,7 @@ MODULE moduleProbe deltaR = NORM2(self%r - part%r) !Only include particle if it is inside the maximum radius - IF (deltaR < self%maxR) THEN + ! IF (deltaR < self%maxR) THEN !find lower index for all dimensions CALL self%findLowerIndex(part%v, i, j, k, inside) @@ -162,40 +162,40 @@ MODULE moduleProbe fk = self%vk(k+1) - part%v(3) fk1 = part%v(3) - self%vk(k) - ! weight = part%weight * DEXP(deltaR/self%maxR) - weight = part%weight + weight = part%weight * DEXP(-deltaR/self%maxR) + ! weight = part%weight !Lock the probe CALL OMP_SET_LOCK(self%lock) !Assign particle weight to distribution function - self%f(i , j , k ) = fi * fj * fk * weight - self%f(i+1, j , k ) = fi1 * fj * fk * weight - self%f(i , j+1, k ) = fi * fj1 * fk * weight - self%f(i+1, j+1, k ) = fi1 * fj1 * fk * weight - self%f(i , j , k+1) = fi * fj * fk1 * weight - self%f(i+1, j , k+1) = fi1 * fj * fk1 * weight - self%f(i , j+1, k+1) = fi * fj1 * fk1 * weight - self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight + self%f(i , j , k ) = self%f(i , j , k ) + fi * fj * fk * weight + self%f(i+1, j , k ) = self%f(i+1, j , k ) + fi1 * fj * fk * weight + self%f(i , j+1, k ) = self%f(i , j+1, k ) + fi * fj1 * fk * weight + self%f(i+1, j+1, k ) = self%f(i+1, j+1, k ) + fi1 * fj1 * fk * weight + self%f(i , j , k+1) = self%f(i , j , k+1) + fi * fj * fk1 * weight + self%f(i+1, j , k+1) = self%f(i+1, j , k+1) + fi1 * fj * fk1 * weight + self%f(i , j+1, k+1) = self%f(i , j+1, k+1) + fi * fj1 * fk1 * weight + self%f(i+1, j+1, k+1) = self%f(i+1, j+1, k+1) + fi1 * fj1 * fk1 * weight !Unlock the probe CALL OMP_UNSET_LOCK(self%lock) END IF - END IF + ! END IF END IF END SUBROUTINE calculate - SUBROUTINE output(self, t) + SUBROUTINE output(self) USE moduleOutput USE moduleRefParam + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(probeDistFunc), INTENT(inout):: self - INTEGER, INTENT(in):: t CHARACTER (LEN=iterationDigits):: tstring CHARACTER (LEN=3):: pstring CHARACTER(:), ALLOCATABLE:: filename @@ -204,14 +204,14 @@ MODULE moduleProbe !Divide by the velocity cube volume self%f = self%f * self%dvInv - WRITE(tstring, iterationFormat) t + WRITE(tstring, iterationFormat) timeStep WRITE(pstring, "(I3.3)") self%id fileName='Probe_' // 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, "(A6, ES15.6E3, A2)") "# t = ", REAL(timeStep)*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) @@ -252,15 +252,15 @@ MODULE moduleProbe END SUBROUTINE doProbes - SUBROUTINE outputProbes(t) + SUBROUTINE outputProbes() + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes IF (probe(i)%update) THEN - CALL probe(i)%output(t) + CALL probe(i)%output() END IF @@ -268,15 +268,15 @@ MODULE moduleProbe END SUBROUTINE outputProbes - SUBROUTINE resetProbes(t) + SUBROUTINE resetProbes() + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes probe(i)%f = 0.D0 - probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0 + probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0 END DO diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index 18fbb7f..e6dc91f 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -160,12 +160,12 @@ MODULE moduleOutput END SUBROUTINE calculateOutput - SUBROUTINE printTime(t, first) + SUBROUTINE printTime(first) USE moduleSpecies USE moduleCompTime + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t LOGICAL, INTENT(in), OPTIONAL:: first CHARACTER(:), ALLOCATABLE:: fileName @@ -187,7 +187,7 @@ MODULE moduleOutput OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField + WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField CLOSE(20) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index a8a5323..740fd14 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -1,52 +1,200 @@ !Module to solve the electromagnetic field MODULE moduleEM + USE moduleMesh + USE moduleTable IMPLICIT NONE - TYPE:: boundaryEM - CHARACTER(:), ALLOCATABLE:: typeEM - INTEGER:: physicalSurface + ! Generic type for electromagnetic boundary conditions + TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric + INTEGER:: nNodes + TYPE(meshNodePointer), ALLOCATABLE:: nodes(:) + + CONTAINS + PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply + !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better. + + END TYPE boundaryEMGeneric + + ABSTRACT INTERFACE + ! Apply boundary condition to the load vector for the Poission equation + SUBROUTINE applyEM_interface(self, vectorF) + IMPORT boundaryEMGeneric + CLASS(boundaryEMGeneric), INTENT(in):: self + REAL(8), INTENT(inout):: vectorF(:) + + END SUBROUTINE applyEM_interface + + END INTERFACE + + TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet REAL(8):: potential CONTAINS - PROCEDURE, PASS:: apply + ! boundaryEMGeneric DEFERRED PROCEDURES + PROCEDURE, PASS:: apply => applyDirichlet - END TYPE boundaryEM + END TYPE boundaryEMDirichlet + + TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime + REAL(8):: potential + TYPE(table1D):: temporalProfile + + CONTAINS + ! boundaryEMGeneric DEFERRED PROCEDURES + PROCEDURE, PASS:: apply => applyDirichletTime + + END TYPE boundaryEMDirichletTime + + ! Container for boundary conditions + TYPE:: boundaryEMCont + CLASS(boundaryEMGeneric), ALLOCATABLE:: obj + + END TYPE boundaryEMCont INTEGER:: nBoundaryEM - TYPE(boundaryEM), ALLOCATABLE:: boundEM(:) + TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:) !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) CONTAINS - !Apply boundary conditions to the K matrix for Poisson's equation - SUBROUTINE apply(self, edge) + SUBROUTINE findNodes(self, physicalSurface) USE moduleMesh IMPLICIT NONE - CLASS(boundaryEM), INTENT(in):: self - CLASS(meshEdge):: edge - INTEGER:: nNodes - INTEGER, ALLOCATABLE:: nodes(:) - INTEGER:: n + CLASS(boundaryEMGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: physicalSurface + CLASS(meshEdge), POINTER:: edge + INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) + INTEGER:: nNodes, nodesNew + INTEGER:: e, n - nNodes = edge%nNodes - nodes = edge%getNodes(nNodes) + !Temporal array to hold nodes + ALLOCATE(nodes(0)) - DO n = 1, nNodes - SELECT CASE(self%typeEM) - CASE ("dirichlet") - mesh%K(nodes(n), :) = 0.D0 - mesh%K(nodes(n), nodes(n)) = 1.D0 - - mesh%nodes(nodes(n))%obj%emData%type = self%typeEM - mesh%nodes(nodes(n))%obj%emData%phi = self%potential + ! Loop thorugh the edges and identify those that are part of the boundary + DO e = 1, mesh%numEdges + edge => mesh%edges(e)%obj + IF (edge%physicalSurface == physicalSurface) THEN + ! Edge is of the right boundary index + ! Get nodes in the edge + nNodes = edge%nNodes + nodesEdge = edge%getNodes(nNodes) + ! Collect all nodes that are not already in the temporal array + DO n = 1, nNodes + IF (ANY(nodes == nodesEdge(n))) THEN + ! Node already in array, skip + CYCLE - END SELECT + ELSE + ! If not, add element to array of nodes + nodes = [nodes, nodesEdge(n)] + + END IF + + END DO + + END IF END DO - END SUBROUTINE apply + ! Point boundary to nodes + nNodes = SIZE(nodes) + ALLOCATE(self%nodes(nNodes)) + self%nNodes = nNodes + DO n = 1, nNodes + self%nodes(n)%obj => mesh%nodes(nodes(n))%obj + + END DO + + END SUBROUTINE findNodes + + ! Initialize Dirichlet boundary condition + SUBROUTINE initDirichlet(self, physicalSurface, potential) + USE moduleRefParam, ONLY: Volt_ref + IMPLICIT NONE + + CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self + INTEGER, INTENT(in):: physicalSurface + REAL(8), INTENT(in):: potential + + ! Allocate boundary edge + ALLOCATE(boundaryEMDirichlet:: self) + + SELECT TYPE(self) + TYPE IS(boundaryEMDirichlet) + self%potential = potential / Volt_ref + + CALL findNodes(self, physicalSurface) + + END SELECT + + END SUBROUTINE initDirichlet + + ! Initialize Dirichlet boundary condition + SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile) + USE moduleRefParam, ONLY: Volt_ref, ti_ref + IMPLICIT NONE + + CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self + INTEGER, INTENT(in):: physicalSurface + REAL(8), INTENT(in):: potential + CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile + + ! Allocate boundary edge + ALLOCATE(boundaryEMDirichletTime:: self) + + SELECT TYPE(self) + TYPE IS(boundaryEMDirichletTime) + self%potential = potential / Volt_ref + + CALL findNodes(self, physicalSurface) + + CALL self%temporalProfile%init(temporalProfile) + + CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0) + + END SELECT + + END SUBROUTINE initDirichletTime + + !Apply Dirichlet boundary condition to the poisson equation + SUBROUTINE applyDirichlet(self, vectorF) + USE moduleMesh + IMPLICIT NONE + + CLASS(boundaryEMDirichlet), INTENT(in):: self + REAL(8), INTENT(inout):: vectorF(:) + INTEGER:: n, ni + + DO n = 1, self%nNodes + self%nodes(n)%obj%emData%phi = self%potential + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi + + END DO + + END SUBROUTINE applyDirichlet + + !Apply Dirichlet boundary condition with time temporal profile + SUBROUTINE applyDirichletTime(self, vectorF) + USE moduleMesh + USE moduleCaseParam, ONLY: timeStep, tauMin + IMPLICIT NONE + + CLASS(boundaryEMDirichletTime), INTENT(in):: self + REAL(8), INTENT(inout):: vectorF(:) + REAL(8):: timeFactor + INTEGER:: n, ni + + timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) + + DO n = 1, self%nNodes + self%nodes(n)%obj%emData%phi = self%potential * timeFactor + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi + + END DO + + END SUBROUTINE applyDirichletTime !Assemble the source vector based on the charge density to solve Poisson's equation SUBROUTINE assembleSourceVector(vectorF, n_e) @@ -60,7 +208,7 @@ MODULE moduleEM REAL(8), ALLOCATABLE:: rho(:) REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes) INTEGER:: nNodes - INTEGER:: e, i, ni + INTEGER:: e, i, ni, b CLASS(meshNode), POINTER:: node ! !$OMP SINGLE diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index c7a0785..f85d812 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -493,52 +493,46 @@ MODULE moduleSolver END SUBROUTINE updateParticleCell !Update the information about if a species needs to be moved this iteration - SUBROUTINE updatePushSpecies(self, t) + SUBROUTINE updatePushSpecies(self) USE moduleSpecies + USE moduleCaseparam, ONLY: timeStep IMPLICIT NONE CLASS(solverGeneric), INTENT(inout):: self - INTEGER, INTENT(in):: t INTEGER:: s DO s=1, nSpecies - self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0 + self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0 END DO END SUBROUTINE updatePushSpecies !Output the different data and information - SUBROUTINE doOutput(t) + SUBROUTINE doOutput() USE moduleMesh USE moduleOutput USE moduleSpecies USE moduleCompTime USE moduleProbe + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t - - IF (t == tInitial) THEN - CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpack_commit.txt') - - END IF - - CALL outputProbes(t) + CALL outputProbes() counterOutput = counterOutput + 1 IF (counterOutput >= triggerOutput .OR. & - t == tFinal .OR. t == tInitial) THEN + timeStep == tFinal .OR. timeStep == tInitial) THEN !Resets output counter counterOutput=0 - CALL mesh%printOutput(t) - IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) - CALL mesh%printEM(t) - WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal + CALL mesh%printOutput() + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl() + CALL mesh%printEM() + WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld - IF (t == 0) THEN + IF (timeStep == 0) THEN WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms" ELSE @@ -556,34 +550,32 @@ MODULE moduleSolver counterCPUTime = counterCPUTime + 1 IF (counterCPUTime >= triggerCPUTime .OR. & - t == tFinal .OR. t == tInitial) THEN + timeStep == tFinal .OR. timeStep == tInitial) THEN !Reset CPU Time counter counterCPUTime = 0 - CALL printTime(t, t == 0) + CALL printTime(timeStep == 0) END IF !Output average values - IF (useAverage .AND. t == tFinal) THEN + IF (useAverage .AND. timeStep == tFinal) THEN CALL mesh%printAverage() END IF END SUBROUTINE doOutput - SUBROUTINE doAverage(t) + SUBROUTINE doAverage() USE moduleAverage USE moduleMesh IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: tAverage, n - IF (useAverage) THEN - tAverage = t - tAverageStart + tAverage = timeStep - tAverageStart IF (tAverage == 1) THEN !First iteration in which average scheme is used