From 9c8bd9a6a5770773696972c4d32d38aa07952a31 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 6 Nov 2007 21:12:04 +0000 Subject: [PATCH] *** empty log message *** --- apps/bvp/AxiStagnBVP.cpp | 38 +- apps/bvp/AxiStagnBVP.h | 5 +- apps/bvp/blasius.mak | 86 +++ apps/bvp/stagnation.cpp | 1143 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 1266 insertions(+), 6 deletions(-) create mode 100644 apps/bvp/blasius.mak create mode 100644 apps/bvp/stagnation.cpp diff --git a/apps/bvp/AxiStagnBVP.cpp b/apps/bvp/AxiStagnBVP.cpp index 36c8f4f16..b094e35f3 100644 --- a/apps/bvp/AxiStagnBVP.cpp +++ b/apps/bvp/AxiStagnBVP.cpp @@ -55,16 +55,18 @@ AxiStagnBVP::AxiStagnBVP(int nsp, int np, double L) : // destructor AxiStagnBVP::~AxiStagnBVP() {} + + // specify guesses for the initial values. These can be anything // that leads to a converged solution. -doublereal initialValue(int n, int j) { +doublereal AxiStagnBVP::initialValue(int n, int j) { switch (n) { case 0: - return m_u0; + return m_uin; case 1: - return m_u0/m_L; + return m_uin/m_L; case 2: - return m_Tinf; + return m_Tin; case 4: return 1.0; default: @@ -72,6 +74,34 @@ doublereal initialValue(int n, int j) { } } +/** + * Set the gas object state to be consistent with the solution at + * point j. + */ +void AxiStagnBVP::setGas(const doublereal* x,int j) { + m_thermo->setTemperature(T(x,j)); + const doublereal* yy = x + m_nv*j + 4; + m_thermo->setMassFractions_NoNorm(yy); + m_thermo->setPressure(m_press); +} + + +/** + * Set the gas state to be consistent with the solution at the + * midpoint between j and j + 1. + */ +void StFlow::setGasAtMidpoint(const doublereal* x,int j) { + m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); + const doublereal* yyj = x + m_nv*j + 4; + const doublereal* yyjp = x + m_nv*(j+1) + 4; + for (int k = 0; k < m_nsp; k++) + m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); + m_thermo->setMassFractions_NoNorm(DATA_PTR(m_ybar)); + m_thermo->setPressure(m_press); +} + + + // Specify the residual. This is where the ODE system and boundary // conditions are specified. The solver will attempt to find a solution // x so that this function returns 0 for all n and j. diff --git a/apps/bvp/AxiStagnBVP.h b/apps/bvp/AxiStagnBVP.h index fc69f4c5a..4f53eb670 100644 --- a/apps/bvp/AxiStagnBVP.h +++ b/apps/bvp/AxiStagnBVP.h @@ -62,12 +62,13 @@ public: // destructor virtual ~AxiStagnBVP() {} + // specify guesses for the initial values. These can be anything // that leads to a converged solution. - virtual doublereal initialValue(int n, int j) { + doublereal AxiStagnBVP::initialValue(int n, int j) { switch (n) { case 0: - return 0.1*z(j); + return m; case 1: return 0.5*z(j); default: diff --git a/apps/bvp/blasius.mak b/apps/bvp/blasius.mak new file mode 100644 index 000000000..310215abe --- /dev/null +++ b/apps/bvp/blasius.mak @@ -0,0 +1,86 @@ +#!/bin/sh + +# This Makefile builds a C++ application that uses Cantera. By +# default, the main program file is 'demo.cpp,' which prints out some +# properties of a reacting gas mixture. + +# To build program 'demo', simply type 'make', or 'make -f ' if this file is named something other than 'Makefile.' + +# Once you have verified that the demo runs, edit this file to replace +# object file 'demo.o' with your own object file or files. + + +#------------------------ edit this block --------------------------------- + +# the name of the executable program to be created +PROG_NAME = blasius.x + +# the object files to be linked together. +OBJS = blasius.o + +# additional flags to be passed to the linker. If your program +# requires other external libraries, put them here +LINK_OPTIONS = -fPIC -L/usr/local/lib -framework Accelerate + +#--------------------------------------------------------------------------- +# You probably don't need to edit anything below. + +# the C++ compiler +CXX = g++ + +# C++ compile flags +CXX_FLAGS = -O3 -Wall -fPIC + +# external libraries +EXT_LIBS = -luser -loneD -lzeroD -lequil -lkinetics -ltransport -lthermo -lctnumerics -lcvode -lctbase -lctmath -ltpx -lctf2c -lconverters -lctcxx + +# Ending C++ linking libraries +LCXX_END_LIBS = -lm + +# the directory where the Cantera libraries are located +CANTERA_LIBDIR=/Applications/Cantera/lib + +# the directory where Cantera include files may be found. +CANTERA_INCDIR=/Applications/Cantera/include + +# flags passed to the C++ compiler/linker for the linking step +LCXXFLAGS = -L$(CANTERA_LIBDIR) -O3 -Wall -fPIC + +# how to compile C++ source files to object files +.cpp.o: + $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + +PROGRAM = $(PROG_NAME)$(EXE_EXT) + +DEPENDS = $(OBJS:.o=.d) + +all: $(PROGRAM) + +$(PROGRAM): $(OBJS) + $(CXX) -o $(PROGRAM) $(OBJS) $(LCXXFLAGS)\ + $(CANTERA_LIBS) $(LINK_OPTIONS) $(EXT_LIBS) \ + $(LCXX_END_LIBS) + +%.d: + g++ -MM -I$(CANTERA_INCDIR) $*.cpp > $*.d + +clean: + $(RM) $(OBJS) $(PROGRAM) + +depends: $(DEPENDS) + cat *.d > .depends + $(RM) $(DEPENDS) + +TAGS: + etags *.h *.cpp + +ifeq ($(wildcard .depends), .depends) +include .depends +endif + + + + + + diff --git a/apps/bvp/stagnation.cpp b/apps/bvp/stagnation.cpp new file mode 100644 index 000000000..e0357b63c --- /dev/null +++ b/apps/bvp/stagnation.cpp @@ -0,0 +1,1143 @@ +/** + * @file AxiStagnBVP.cpp + */ + +/* + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2002 California Institute of Technology + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#pragma warning(disable:4267) +#endif + +#include +#include + +#include "AxiStagnBVP.h" +#include "ArrayViewer.h" +#include "ctml.h" +#include "MultiJac.h" + +using namespace ctml; +using namespace Cantera; +using namespace std; + +static void st_drawline() { + writelog("\n-------------------------------------" + "------------------------------------------"); +} + +AxiStagnBVP::AxiStagnBVP(igthermo_t* ph, int nsp, int points) : + Domain1D(nsp+4, points), + m_inlet_u(0.0), + m_inlet_V(0.0), + m_inlet_T(-1.0), + m_surface_T(-1.0), + m_press(-1.0), + m_nsp(nsp), + m_thermo(0), + m_kin(0), + m_trans(0), + m_jac(0), + m_ok(false), + m_do_soret(false), + m_transport_option(-1), + m_efctr(0.0) + { + m_type = cFlowType; + + m_points = points; + m_thermo = ph; + + if (ph == 0) return; // used to create a dummy object + + int nsp2 = m_thermo->nSpecies(); + if (nsp2 != m_nsp) { + m_nsp = nsp2; + Domain1D::resize(m_nsp+4, points); + } + + + // make a local copy of the species molecular weight vector + m_wt = m_thermo->molecularWeights(); + + // the species mass fractions are the last components in the solution + // vector, so the total number of components is the number of species + // plus the offset of the first mass fraction. + m_nv = c_offset_Y + m_nsp; + + // enable all species equations by default + m_do_species.resize(m_nsp, true); + + // but turn off the energy equation at all points + m_do_energy.resize(m_points,false); + + m_diff.resize(m_nsp*m_points); + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_surfdot.resize(m_nsp, 0.0); + m_ybar.resize(m_nsp); + + + //-------------- default solution bounds -------------------- + + vector_fp vmin(m_nv), vmax(m_nv); + + // no bounds on u + vmin[0] = -1.e20; + vmax[0] = 1.e20; + + // V + vmin[1] = -1.e20; + vmax[1] = 1.e20; + + // temperature bounds + vmin[2] = 200.0; + vmax[2]= 1.e9; + + // lamda should be negative + vmin[3] = -1.e20; + vmax[3] = 1.e20; + + // mass fraction bounds + int k; + for (k = 0; k < m_nsp; k++) { + vmin[4+k] = -1.0e-5; + vmax[4+k] = 1.0e5; + } + setBounds(vmin.size(), DATA_PTR(vmin), vmax.size(), DATA_PTR(vmax)); + + + //-------------------- default error tolerances ---------------- + vector_fp rtol(m_nv, 1.0e-8); + vector_fp atol(m_nv, 1.0e-15); + setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),false); + setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),true); + + //-------------------- grid refinement ------------------------- + m_refiner->setActive(0, false); + m_refiner->setActive(1, false); + m_refiner->setActive(2, false); + m_refiner->setActive(3, false); + + vector_fp gr; + for (int ng = 0; ng < m_points; ng++) gr.push_back(1.0*ng/m_points); + setupGrid(m_points, DATA_PTR(gr)); + setID("stagnation flow"); + } + + + /** + * Change the grid size. Called after grid refinement. + */ + void AxiStagnBVP::resize(int ncomponents, int points) { + Domain1D::resize(ncomponents, points); + m_rho.resize(m_points, 0.0); + m_wtm.resize(m_points, 0.0); + m_cp.resize(m_points, 0.0); + m_enth.resize(m_points, 0.0); + m_visc.resize(m_points, 0.0); + m_tcon.resize(m_points, 0.0); + + if (m_transport_option == c_Mixav_Transport) { + m_diff.resize(m_nsp*m_points); + } + else { + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_diff.resize(m_nsp*m_points); + } + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_do_energy.resize(m_points,false); + + m_fixedy.resize(m_nsp, m_points); + m_fixedtemp.resize(m_points); + + m_dz.resize(m_points-1); + m_z.resize(m_points); + } + + + void AxiStagnBVP::setupGrid(int n, const doublereal* z) { + resize(m_nv, n); + int j; + + m_z[0] = z[0]; + for (j = 1; j < m_points; j++) { + m_z[j] = z[j]; + m_dz[j-1] = m_z[j] - m_z[j-1]; + } + } + + + /** + * Install a transport manager. + */ + void AxiStagnBVP::setTransport(Transport& trans, bool withSoret) { + m_trans = &trans; + m_do_soret = withSoret; + + if (m_trans->model() == cMulticomponent) { + m_transport_option = c_Multi_Transport; + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_diff.resize(m_nsp*m_points); + m_dthermal.resize(m_nsp, m_points, 0.0); + } + else if (m_trans->model() == cMixtureAveraged) { + m_transport_option = c_Mixav_Transport; + m_diff.resize(m_nsp*m_points); + if (withSoret) + throw CanteraError("setTransport", + "Thermal diffusion (the Soret effect) " + "requires using a multicomponent transport model."); + } + else + throw CanteraError("setTransport","unknown transport model."); + } + + void AxiStagnBVP::enableSoret(bool withSoret) { + if (m_transport_option == c_Multi_Transport) + m_do_soret = withSoret; + else { + throw CanteraError("setTransport", + "Thermal diffusion (the Soret effect) " + "requires using a multicomponent transport model."); + } + } + + + /** + * Set the gas object state to be consistent with the solution at + * point j. + */ + void AxiStagnBVP::setGas(const doublereal* x,int j) { + m_thermo->setTemperature(T(x,j)); + const doublereal* yy = x + m_nv*j + c_offset_Y; + m_thermo->setMassFractions_NoNorm(yy); + m_thermo->setPressure(m_press); + } + + + /** + * Set the gas state to be consistent with the solution at the + * midpoint between j and j + 1. + */ + void AxiStagnBVP::setGasAtMidpoint(const doublereal* x,int j) { + m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); + const doublereal* yyj = x + m_nv*j + c_offset_Y; + const doublereal* yyjp = x + m_nv*(j+1) + c_offset_Y; + for (int k = 0; k < m_nsp; k++) + m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); + m_thermo->setMassFractions_NoNorm(DATA_PTR(m_ybar)); + m_thermo->setPressure(m_press); + } + + + void AxiStagnBVP::_finalize(const doublereal* x) { + int k, j; + doublereal zz, tt; + int nz = m_zfix.size(); + bool e = m_do_energy[0]; + for (j = 0; j < m_points; j++) { + if (e || nz == 0) + setTemperature(j, T(x, j)); + else { + zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); + tt = linearInterp(zz, m_zfix, m_tfix); + setTemperature(j, tt); + } + for (k = 0; k < m_nsp; k++) { + setMassFraction(j, k, Y(x, k, j)); + } + } + if (e) solveEnergyEqn(); + } + + + //------------------------------------------------------ + + /** + * Evaluate the residual function for axisymmetric stagnation + * flow. If jpt is less than zero, the residual function is + * evaluated at all grid points. If jpt >= 0, then the residual + * function is only evaluated at grid points jpt-1, jpt, and + * jpt+1. This option is used to efficiently evaluate the + * Jacobian numerically. + * + */ + + void AxiStagnFlowBVP::prepare(doublereal* x) { + int j; + + // update thermo properties + for (j = 0; j < m_points; j++) { + setGas(x,j); + m_rho[j] = m_thermo->density(); + m_wtm[j] = m_thermo->meanMolecularWeight(); + m_cp[j] = m_thermo->cp_mass(); + } + + // update transport properties + int k,m; + + if (m_transport_option == c_Mixav_Transport) { + for (j = 0; j < m_points; j++) { + setGasAtMidpoint(x,j); + m_visc[j] = m_trans->viscosity(); + m_trans->getMixDiffCoeffs(DATA_PTR(m_diff) + j*m_nsp); + m_tcon[j] = m_trans->thermalConductivity(); + } + } + else if (m_transport_option == c_Multi_Transport) { + doublereal sum, sumx, wtm, dz; + doublereal eps = 1.0e-12; + for (m = 0; m < m_points-1; m++) { + setGasAtMidpoint(x,m); + dz = m_z[m+1] - m_z[m]; + wtm = m_thermo->meanMolecularWeight(); + + m_visc[m] = m_trans->viscosity(); + m_trans->getMultiDiffCoeffs(m_nsp, + DATA_PTR(m_multidiff) + mindex(0,0,m)); + + for (k = 0; k < m_nsp; k++) { + sum = 0.0; + sumx = 0.0; + for (j = 0; j < m_nsp; j++) { + if (j != k) { + sum += m_wt[j]*m_multidiff[mindex(k,j,m)]* + ((X(x,j,m+1) - X(x,j,m))/dz + eps); + sumx += (X(x,j,m+1) - X(x,j,m))/dz; + } + } + m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps)); + } + + m_tcon[m] = m_trans->thermalConductivity(); + if (m_do_soret) { + m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp); + } + } + } + } + + + void AxiStagnFlowBVP::residual(doublereal* x, + int n, int j) { + + int j, k; + + updateDiffFluxes(x, j0, j1); + + + //---------------------------------------------------- + // evaluate the residual equations at all required + // grid points + //---------------------------------------------------- + + doublereal sum, sum2, dtdzj; + + + //---------------------------------------------- + // left boundary + //---------------------------------------------- + + if (j == 0) { + + // Continuity. This propagates information right-to-left, + // since rho_u at point 0 is dependent on rho_u at point 1, + // but not on mdot from the inlet. + if (n == c_offset_U) { + return -(rho_u(x,1) - rho_u(x,0))/m_dz[0] - (density(1)*V(x,1) + density(0)*V(x,0)); + } + + // the inlet (or other) object connected to this one + // will modify these equations by subtracting its values + // for V, T, and mdot. As a result, these residual equations + // will force the solution variables to the values for + // the boundary object + else if (n == c_offset_V) return V(x,0); + else if (n == c_offset_T) return m_Tsurf - T(x,0); + rsd[index(c_offset_T,0)] = T(x,0); + rsd[index(c_offset_L,0)] = -rho_u(x,0); + + // The default boundary condition for species is zero + // flux. However, the boundary object may modify + // this. + sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,0); + rsd[index(c_offset_Y + k, 0)] = + -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); + } + rsd[index(c_offset_Y, 0)] = 1.0 - sum; + } + + + //---------------------------------------------- + // + // right boundary + // + //---------------------------------------------- + + else if (j == m_points - 1) { + + // the boundary object connected to the right of this + // one may modify or replace these equations. The + // default boundary conditions are zero u, V, and T, + // and zero diffusive flux for all species. + + rsd[index(0,j)] = rho_u(x,j); + rsd[index(1,j)] = V(x,j); + rsd[index(2,j)] = T(x,j); + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + doublereal sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); + } + rsd[index(4,j)] = 1.0 - sum; + diag[index(4,j)] = 0; + + } + + + //------------------------------------------ + // interior points + //------------------------------------------ + + else { + + //---------------------------------------------- + // Continuity equation + // + // Note that this propagates the mass flow rate + // information to the left (j+1 -> j) from the + // value specified at the right boundary. The + // lambda information propagates in the opposite + // direction. + // + // d(\rho u)/dz + 2\rho V = 0 + // + //------------------------------------------------ + + rsd[index(c_offset_U,j)] = + -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); + + //algebraic constraint + diag[index(c_offset_U, j)] = 0; + + + //------------------------------------------------ + // Radial momentum equation + // + // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda + // + //------------------------------------------------- + rsd[index(c_offset_V,j)] + = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) + - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] + - rdt*(V(x,j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; + + + //------------------------------------------------- + // Species equations + // + // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k + // + //------------------------------------------------- + getWdot(x,j); + + doublereal convec, diffus; + for (k = 0; k < m_nsp; k++) { + convec = rho_u(x,j)*dYdz(x,k,j); + diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + /(z(j+1) - z(j-1)); + rsd[index(c_offset_Y + k, j)] + = (m_wt[k]*(wdot(k,j) ) + - convec - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; + } + + + //----------------------------------------------- + // energy equation + //----------------------------------------------- + + if (m_do_energy[j]) { + + setGas(x,j); + + // heat release term + const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); + const vector_fp& cp_R = m_thermo->cp_R_ref(); + + sum = 0.0; + sum2 = 0.0; + doublereal flxk; + for (k = 0; k < m_nsp; k++) { + flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + sum += wdot(k,j)*h_RT[k]; + sum2 += flxk*cp_R[k]/m_wt[k]; + } + sum *= GasConstant * T(x,j); + dtdzj = dTdz(x,j); + sum2 *= GasConstant * dtdzj; + + rsd[index(c_offset_T, j)] = + - m_cp[j]*rho_u(x,j)*dtdzj + - divHeatFlux(x,j) - sum - sum2; + rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); + + rsd[index(c_offset_T, j)] = + rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j)); + + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + diag[index(c_offset_T, j)] = 1; + } + + // residual equations if the energy equation is disabled + + if (!m_do_energy[j]) { + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); + diag[index(c_offset_T, j)] = 0; + } + + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + } + } + } + + + + /** + * Update the transport properties at grid points in the range + * from j0 to j1, based on solution x. + */ + void AxiStagnBVP::updateTransport(doublereal* x,int j0, int j1) { + int j,k,m; + + if (m_transport_option == c_Mixav_Transport) { + for (j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); + m_trans->getMixDiffCoeffs(DATA_PTR(m_diff) + j*m_nsp); + m_tcon[j] = m_trans->thermalConductivity(); + } + } + else if (m_transport_option == c_Multi_Transport) { + doublereal sum, sumx, wtm, dz; + doublereal eps = 1.0e-12; + for (m = j0; m < j1; m++) { + setGasAtMidpoint(x,m); + dz = m_z[m+1] - m_z[m]; + wtm = m_thermo->meanMolecularWeight(); + + m_visc[m] = (m_dovisc ? m_trans->viscosity() : 0.0); + + m_trans->getMultiDiffCoeffs(m_nsp, + DATA_PTR(m_multidiff) + mindex(0,0,m)); + + for (k = 0; k < m_nsp; k++) { + sum = 0.0; + sumx = 0.0; + for (j = 0; j < m_nsp; j++) { + if (j != k) { + sum += m_wt[j]*m_multidiff[mindex(k,j,m)]* + ((X(x,j,m+1) - X(x,j,m))/dz + eps); + sumx += (X(x,j,m+1) - X(x,j,m))/dz; + } + } + m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps)); + } + + m_tcon[m] = m_trans->thermalConductivity(); + if (m_do_soret) { + m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp); + } + } + } + } + + + //------------------------------------------------------ + + /** + * Evaluate the residual function for axisymmetric stagnation + * flow. If jpt is less than zero, the residual function is + * evaluated at all grid points. If jpt >= 0, then the residual + * function is only evaluated at grid points jpt-1, jpt, and + * jpt+1. This option is used to efficiently evaluate the + * Jacobian numerically. + * + */ + + void FreeFlame::eval(int jg, doublereal* xg, + doublereal* rg, integer* diagg, doublereal rdt) { + + // if evaluating a Jacobian, and the global point is outside + // the domain of influence for this domain, then skip + // evaluating the residual + if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; + + // if evaluating a Jacobian, compute the steady-state residual + if (jg >= 0) rdt = 0.0; + + // start of local part of global arrays + doublereal* x = xg + loc(); + doublereal* rsd = rg + loc(); + integer* diag = diagg + loc(); + + int jmin, jmax, jpt; + jpt = jg - firstPoint(); + + if (jg < 0) { // evaluate all points + jmin = 0; + jmax = m_points - 1; + } + else { // evaluate points for Jacobian + jmin = max(jpt-1, 0); + jmax = min(jpt+1,m_points-1); + } + + // properties are computed for grid points from j0 to j1 + int j0 = max(jmin-1,0); + int j1 = min(jmax+1,m_points-1); + + + int j, k; + + + //----------------------------------------------------- + // update properties + //----------------------------------------------------- + + // update thermodynamic properties only if a Jacobian is not + // being evaluated + if (jpt < 0) { + updateThermo(x, j0, j1); + updateTransport(x, j0, j1); + } + + // update the species diffusive mass fluxes whether or not a + // Jacobian is being evaluated + updateDiffFluxes(x, j0, j1); + + + //---------------------------------------------------- + // evaluate the residual equations at all required + // grid points + //---------------------------------------------------- + + doublereal sum, sum2, dtdzj; + + for (j = jmin; j <= jmax; j++) { + + + //---------------------------------------------- + // left boundary + //---------------------------------------------- + + if (j == 0) { + + // these may be modified by a boundary object + + // Continuity. This propagates information right-to-left, + // since rho_u at point 0 is dependent on rho_u at point 1, + // but not on mdot from the inlet. + rsd[index(c_offset_U,0)] = + -(rho_u(x,1) - rho_u(x,0))/m_dz[0] + -(density(1)*V(x,1) + density(0)*V(x,0)); + + // the inlet (or other) object connected to this one + // will modify these equations by subtracting its values + // for V, T, and mdot. As a result, these residual equations + // will force the solution variables to the values for + // the boundary object + rsd[index(c_offset_V,0)] = V(x,0); + rsd[index(c_offset_T,0)] = T(x,0); + rsd[index(c_offset_L,0)] = -rho_u(x,0); + + // The default boundary condition for species is zero + // flux + sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,0); + rsd[index(c_offset_Y + k, 0)] = + -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); + } + rsd[index(c_offset_Y, 0)] = 1.0 - sum; + } + + + //---------------------------------------------- + // + // right boundary + // + //---------------------------------------------- + + else if (j == m_points - 1) { + + // the boundary object connected to the right of this + // one may modify or replace these equations. The + // default boundary conditions are zero u, V, and T, + // and zero diffusive flux for all species. + + // zero gradient + rsd[index(0,j)] = rho_u(x,j) - rho_u(x,j-1); + rsd[index(1,j)] = V(x,j); + rsd[index(2,j)] = T(x,j) - T(x,j-1); + doublereal sum = 0.0; + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); + } + rsd[index(4,j)] = 1.0 - sum; + diag[index(4,j)] = 0; + } + + //------------------------------------------ + // interior points + //------------------------------------------ + + else { + + //---------------------------------------------- + // Continuity equation + //---------------------------------------------- + + if (grid(j) > m_zfixed){ + rsd[index(c_offset_U,j)] = + - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1] + - (density(j-1)*V(x,j-1) + density(j)*V(x,j)); + } + + else if (grid(j) == m_zfixed){ + if (m_do_energy[j]) { + rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + } + else { + rsd[index(c_offset_U,j)] = (rho_u(x,j) + - m_rho[0]*0.3); + } + } + else if(grid(j) < m_zfixed){ + rsd[index(c_offset_U,j)] = + - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + - (density(j+1)*V(x,j+1) + density(j)*V(x,j)); + } + //algebraic constraint + diag[index(c_offset_U, j)] = 0; + + //------------------------------------------------ + // Radial momentum equation + // + // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda + // + //------------------------------------------------- + rsd[index(c_offset_V,j)] + = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) + - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] + - rdt*(V(x,j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; + + + //------------------------------------------------- + // Species equations + // + // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k + // + //------------------------------------------------- + getWdot(x,j); + + doublereal convec, diffus; + for (k = 0; k < m_nsp; k++) { + convec = rho_u(x,j)*dYdz(x,k,j); + diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + /(z(j+1) - z(j-1)); + rsd[index(c_offset_Y + k, j)] + = (m_wt[k]*(wdot(k,j) ) + - convec - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; + } + + + //----------------------------------------------- + // energy equation + //----------------------------------------------- + + if (m_do_energy[j]) { + + setGas(x,j); + + // heat release term + const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); + const vector_fp& cp_R = m_thermo->cp_R_ref(); + + sum = 0.0; + sum2 = 0.0; + doublereal flxk; + for (k = 0; k < m_nsp; k++) { + flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + sum += wdot(k,j)*h_RT[k]; + sum2 += flxk*cp_R[k]/m_wt[k]; + } + sum *= GasConstant * T(x,j); + dtdzj = dTdz(x,j); + sum2 *= GasConstant * dtdzj; + + rsd[index(c_offset_T, j)] = + - m_cp[j]*rho_u(x,j)*dtdzj + - divHeatFlux(x,j) - sum - sum2; + rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); + + rsd[index(c_offset_T, j)] = + rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j)); + + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + diag[index(c_offset_T, j)] = 1; + } + // residual equations if the energy equation is disabled + else { + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); + diag[index(c_offset_T, j)] = 0; + } + + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + } + } + } + + + /** + * Print the solution. + */ + void AxiStagnBVP::showSolution(const doublereal* x) { + int nn = m_nv/5; + int i, j, n; + //char* buf = new char[100]; + char buf[100]; + + // The mean molecular weight is needed to convert + updateThermo(x, 0, m_points-1); + + sprintf(buf, " Pressure: %10.4g Pa \n", m_press); + writelog(buf); + for (i = 0; i < nn; i++) { + st_drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10s ",componentName(i*5 + n).c_str()); + writelog(buf); + } + st_drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10.4g ",component(x, i*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } + int nrem = m_nv - 5*nn; + st_drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10s ", componentName(nn*5 + n).c_str()); + writelog(buf); + } + st_drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10.4g ",component(x, nn*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } + + + /** + * Update the diffusive mass fluxes. + */ + void AxiStagnBVP::updateDiffFluxes(const doublereal* x, int j0, int j1) { + int j, k, m; + doublereal sum, wtm, rho, dz, gradlogT; + + switch (m_transport_option) { + + case c_Mixav_Transport: + case c_Multi_Transport: + for (j = j0; j < j1; j++) { + sum = 0.0; + wtm = m_wtm[j]; + rho = density(j); + dz = z(j+1) - z(j); + + for (k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); + m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; + sum -= m_flux(k,j); + } + // correction flux to insure that \sum_k Y_k V_k = 0. + for (k = 0; k < m_nsp; k++) m_flux(k,j) += sum*Y(x,k,j); + } + break; + + default: + throw CanteraError("updateDiffFluxes","unknown transport model"); + } + + if (m_do_soret) { + for (m = j0; m < j1; m++) { + gradlogT = 2.0*(T(x,m+1) - T(x,m))/(T(x,m+1) + T(x,m)); + for (k = 0; k < m_nsp; k++) { + m_flux(k,m) -= m_dthermal(k,m)*gradlogT; + } + } + } + } + + + string AxiStagnBVP::componentName(int n) const { + switch(n) { + case 0: return "u"; + case 1: return "V"; + case 2: return "T"; + case 3: return "lambda"; + default: + if (n >= (int) c_offset_Y && n < (int) (c_offset_Y + m_nsp)) { + return m_thermo->speciesName(n - c_offset_Y); + } + else + return ""; + } + } + + + //added by Karl Meredith + int AxiStagnBVP::componentIndex(string name) const { + + + if(name=="u") {return 0;} + else if (name=="V") {return 1;} + else if (name=="T") {return 2;} + else if (name=="lambda") {return 3;} + else { + for (int n=4;n ignored; + int nsp = m_thermo->nSpecies(); + vector_int did_species(nsp, 0); + + vector str; + dom.getChildren("string",str); + int nstr = static_cast(str.size()); + for (int istr = 0; istr < nstr; istr++) { + const XML_Node& nd = *str[istr]; + writelog(nd["title"]+": "+nd.value()+"\n"); + } + + //map params; + double pp = -1.0; + pp = getFloat(dom, "pressure", "pressure"); + setPressure(pp); + + + vector d; + dom.child("grid_data").getChildren("floatArray",d); + int nd = static_cast(d.size()); + + vector_fp x; + int n, np = 0, j, ks, k; + string nm; + bool readgrid = false, wrote_header = false; + for (n = 0; n < nd; n++) { + const XML_Node& fa = *d[n]; + nm = fa["title"]; + if (nm == "z") { + getFloatArray(fa,x,false); + np = x.size(); + writelog("Grid contains "+int2str(np)+ + " points.\n"); + readgrid = true; + setupGrid(np, DATA_PTR(x)); + } + } + if (!readgrid) { + throw CanteraError("AxiStagnBVP::restore", + "domain contains no grid points."); + } + + writelog("Importing datasets:\n"); + for (n = 0; n < nd; n++) { + const XML_Node& fa = *d[n]; + nm = fa["title"]; + getFloatArray(fa,x,false); + if (nm == "u") { + writelog("axial velocity "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) { + soln[index(0,j)] = x[j]; + } + } + else { + goto error; + } + } + else if (nm == "z") { + ; // already read grid + } + else if (nm == "V") { + writelog("radial velocity "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(1,j)] = x[j]; + } + else goto error; + } + else if (nm == "T") { + writelog("temperature "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(2,j)] = x[j]; + + // For fixed-temperature simulations, use the + // imported temperature profile by default. If + // this is not desired, call setFixedTempProfile + // *after* restoring the solution. + + vector_fp zz(np); + for (int jj = 0; jj < np; jj++) + zz[jj] = (grid(jj) - zmin())/(zmax() - zmin()); + setFixedTempProfile(zz, x); + } + else goto error; + } + else if (nm == "L") { + writelog("lambda "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(3,j)] = x[j]; + } + else goto error; + } + else if (m_thermo->speciesIndex(nm) >= 0) { + writelog(nm+" "); + if ((int) x.size() == np) { + k = m_thermo->speciesIndex(nm); + did_species[k] = 1; + for (j = 0; j < np; j++) + soln[index(k+4,j)] = x[j]; + } + } + else + ignored.push_back(nm); + } + + if (ignored.size() != 0) { + writelog("\n\n"); + writelog("Ignoring datasets:\n"); + int nn = static_cast(ignored.size()); + for (int n = 0; n < nn; n++) { + writelog(ignored[n]+" "); + } + } + + for (ks = 0; ks < nsp; ks++) { + if (did_species[ks] == 0) { + if (!wrote_header) { + writelog("Missing data for species:\n"); + wrote_header = true; + } + writelog(m_thermo->speciesName(ks)+" "); + } + } + + return; + error: + throw CanteraError("AxiStagnBVP::restore","Data size error"); + } + + + + void AxiStagnBVP::save(XML_Node& o, doublereal* sol) { + int k; + + ArrayViewer soln(m_nv, m_points, sol + loc()); + + XML_Node& flow = (XML_Node&)o.addChild("domain"); + flow.addAttribute("type",flowType()); + flow.addAttribute("id",m_id); + flow.addAttribute("points",m_points); + flow.addAttribute("components",m_nv); + + if (m_desc != "") addString(flow,"description",m_desc); + XML_Node& gv = flow.addChild("grid_data"); + addFloat(flow, "pressure", m_press, "Pa", "pressure"); + addFloatArray(gv,"z",m_z.size(),DATA_PTR(m_z), + "m","length"); + vector_fp x(static_cast(soln.nColumns())); + + soln.getRow(0,DATA_PTR(x)); + addFloatArray(gv,"u",x.size(),DATA_PTR(x),"m/s","velocity"); + + soln.getRow(1,DATA_PTR(x)); + addFloatArray(gv,"V", + x.size(),DATA_PTR(x),"1/s","rate"); + + soln.getRow(2,DATA_PTR(x)); + addFloatArray(gv,"T",x.size(),DATA_PTR(x),"K","temperature",0.0); + + soln.getRow(3,DATA_PTR(x)); + addFloatArray(gv,"L",x.size(),DATA_PTR(x),"N/m^4"); + + for (k = 0; k < m_nsp; k++) { + soln.getRow(4+k,DATA_PTR(x)); + addFloatArray(gv,m_thermo->speciesName(k), + x.size(),DATA_PTR(x),"","massFraction",0.0,1.0); + } + } + + + void AxiStagnBVP::setJac(MultiJac* jac) { + m_jac = jac; + } + + +} // namespace