guyer
commented
11 months ago -
FiPy supports lazy evaluation and can operate on whole fields at once, so, rather than calculating a field and assigning it to the value of a
CellVariable, you should write the expression you want and us it as aCellVariable. Here's my attempt to do so withcompute_species_source_term():def compute_diffusion_coefficient( Z_st, SDR_st, mesh ): Z = mesh.faceCenters[0] return SDR_st * np.exp(2. * (np.square(erfcinv(2. * Z_st)) - np.square(erfcinv(2. * Z)))) / 2. def compute_species_source_term( phase, phi_vec, mesh, species_index_of_interest ): # update thermophysical state in terms of mass fractions and temperature phase.Y_vec = phi_vec[:len(phase.species_name_vec)] phase.T = phi_vec[-1] return phase.net_mass_production_rate_vec()[species_index_of_interest] / phase.rho() def compute_energy_source_term( phase, phi_vec, mesh ): """ Description: Source term in energy conservation equation. """ # update thermophysical state in terms of mass fractions and temperature phase.Y_vec = phi_vec[:len(phase.species_name_vec)] phase.T = phi_vec[-1] # compute auxiliary [J m^-3 s^-1] vector auxiliary_vec = phase.ha_vec().dot(phase.net_mass_production_rate_vec()) # store respective value of energy source term in each cell center return -1. / (phase.rho() * phase.cp()) * auxiliary_vecSimilar changes are likely needed in
phase.net_mass_production_rate_vec(). -
The documentation does not say that no preconditioners are implemented for PETSc. It says
FiPy does not implement any precoditioner objects for PETSc. Simply pass one of the PCType strings in the precon= argument when declaring the solver.
So, you can, for example, call
LinearGMRESSolver(..., precon="jacobi")orLinearGMRESSolver(..., precon="rowscalingviennacl")to get the corresponding preconditioner from https://petsc.org/release/manualpages/PC/PCType/. (No, I have no idea what the latter one is).
marvinpom
Hello,
I have some issues with applying FiPy to my problem and would like to ask for your assistance.
Physical problem description:
I would like to use FiPy to solve the steady-state counterflow diffusion flame problem at constant prescribed background pressure, i.e. the classical "flamelet" equations in mixture fraction coordinates (1D mesh) for combustion modeling. The coupled set of equations is composed of the species transport equations as well as the energy equation:
$\frac{\partial Y{i}}{\partial t} = 0 = \frac{\chi}{2}\frac{\partial^{2} Y{i}}{\partial Z^{2}} + \frac{\dot{\omega{i}}}{\rho}$ where $i = 1,...,N{\mathrm{species}}$ $\frac{\partial T}{\partial t} = 0 = \frac{\chi}{2}\frac{\partial^{2} T}{\partial Z^{2}} - \frac{1}{\rho c{p}} \sum{i}{h{i} \dot{\omega}{i}}$
$Y{i}$: mass fraction of chemical species $i$ [kg kg^-1] $T$: temperature [K] $\chi$: scalar dissipation rate [s^-1] $Z$: mixture fraction [-] $\dot{\omega{i}}$: chemical source term of chemical species $i$ [kg m^-3 s^-1] $\rho$: mass-specific density [kg m^-3] $c{p}$: mass-specific isobaric heat capacity [J kg^-1 K^-1] $h{i}$: mass-specific (absolute) enthalpy of chemical species $i$ [J kg^-1]
In the end, I am interested in the chemical mixture composition in terms of mass fractions as a function of mixture fraction $Y_{i}(Z)$ as well as in the temperature profile $T(Z)$. The mixture fraction represents the local fuel content in the mixture and is restricted to the interval [0, 1]. The scalar dissipation rate, which acts as the diffusion coefficient of the problem under consideration ($D = \chi/2$), depends on the mesh (i.e. the mixture fraction) and results from the analytical profile
$\chi(Z) = \chi{\mathrm{st}}\mathrm{exp}(2((\mathrm{erfc}^{-1}(2 Z{\mathrm{st}}))^{2}-(\mathrm{erfc}^{-1}(2 Z))^{2}))$
where the index "st" represents the value at stoichiometric mixing conditions. The value of $\chi{\mathrm{st}}$ is prescribed and kept constant. The value of $Z{\mathrm{st}}$ is also constant for a fixed propellant combination and thus known a-priori.
The remaining variables are functions of pressure $p$ and the solution variables $Y{i}, T$, i.e. $\dot{\omega{i}}, \rho, c{p}, h{i} = f(p, Y_{i}, T)$. I can determine their values with the help of another class of my python program, as soon as I know (the pressure,) the mass fraction vector and the temperature.
I would like to apply Dirichlet boundary conditions for the propellant compositions and injection temperatures:
@ oxidizer inlet (corresponds to mixture fraction $Z = 0$ ~ no fuel in mixture) if i == oxindex: $Y{i}(Z=0) = 1$ else: $Y_{i}(Z=0) = 0$
$T(Z=0) = T_{\mathrm{ox}}$
@ fuel inlet (corresponds to mixture fraction $Z = 1$ ~ pure fuel in mixture) if i == fuindex: $Y{i}(Z=1) = 1$ else: $Y_{i}(Z=1) = 0$
$T(Z=1) = T_{\mathrm{fu}}$
Furthermore, resonable initial guesses for the profiles of the solution variables are available.
I would like to employ the Generalized Minimum Residual (GMRES) solver method from PETSc.
Code:
I attach the code snippets of my attempt to solve the Flamelet equations with FiPy as .txt file:
solver.txt
Issues:
With this code, I observe a strange solver behavior, which may be due to the fact that both the respective source terms in the species transport equations $S{\mathrm{species}} = \frac{\dot{\omega{i}}}{\rho}$ as well as the source term in the energy equation $S{\mathrm{energy}} = -\frac{1}{\rho c{p}} \sum{i}{h{i} \dot{\omega}_{i}}$ are evaluated only once at the beginning and are not updated again according to the current values of the solution variables in each iteration. How can I ensure that the source terms always correspond to the current values of the solution variables and are therefore consistent?
The documentation (https://www.ctcms.nist.gov/fipy/documentation/SOLVERS.html) states that no preconditioners are implemented for PETSc. However, in the generated log file it is stated that the incomplete LU-decomposition ("ilu") is applied. Is the documentation not up to date or is it just a "dummy" entry in the log-file? Since the system of equations is very stiff due to the chemical reaction terms, a suitable preconditioning is probably indispensable. Since Trilinos and PySparse are only available for python 2 and the documentation states that there are no preconditioners available for either PETSc or scipy.sparse solvers as well, I wonder what I could do. Is PyAMG in combination with a scipy solver the only alternative?
I would like to thank you for your help in advance.
Best regards