Switch updateElectronTransport(true);

Info<< "Reading physicalProperties\n" << endl;

IOdictionary physicalProperties
(
    IOobject
    (
        "physicalProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

Switch constETransport
(
    physicalProperties.lookupOrDefault("constETransport", false)
);

scalar relaxDrift
(
    physicalProperties.lookupOrDefault("relaxDrift", 1.0)
);

dimensionedScalar epsilon0
(
    physicalProperties.lookup("epsilon0")
);

Switch mobility_f_of_Te = physicalProperties.lookupOrDefault("mobility_f_of_Te", false);

Switch calculateDe = physicalProperties.lookupOrDefault("calculateDe", false);

// Convert E/n in SI unit to table unit.
// Default V m^2 => Td
scalar TeToTableUnit (
        physicalProperties.lookupOrDefault("TeToTableUnit", 1.0)
        );

// Convert E/n in SI unit to table unit.
// Default V m^2 => Td
scalar EnToTableUnit (
        physicalProperties.lookupOrDefault("EnToTableUnit", 1.0e21)
        );


autoPtr< DataEntry< scalar > > pmueN (
        DataEntry<scalar>::New("mueN", physicalProperties));
const DataEntry<scalar> &mueN(pmueN());

// Convert mu_e * n_g value from table into SI unit.
scalar mueNFac (
        physicalProperties.lookupOrDefault("mueNFac", 1.0)
        );


autoPtr< DataEntry< scalar > > pDeN (
        DataEntry<scalar>::New("DeN", physicalProperties));
const DataEntry<scalar> &DeN(pDeN());

// Convert D_e * n_g value from table into SI unit.
scalar DeNFac (
        physicalProperties.lookupOrDefault("DeNFac", 1.0)
        );


autoPtr< DataEntry< scalar > > pTeOfEn(
        DataEntry<scalar>::New ("TeOfEn", physicalProperties));
const DataEntry<scalar> &TeOfEn(pTeOfEn());

// Convert T_e value from table into SI unit.
scalar TeFac (
        physicalProperties.lookupOrDefault("TeFac", 1.0)
        );


Info<< "Reading field Phi\n" << endl;
volScalarField Phi
(
    IOobject
    (
        "Phi",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Creating field electric field\n" << endl;
volVectorField E
(
    IOobject
    (
        "E",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    -fvc::grad(Phi)
);

surfaceScalarField snE
(
    IOobject
    (
        "snE",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    -fvc::snGrad(Phi)
);


Info<< "Creating reaction model\n" << endl;

autoPtr<combustionModels::psiCombustionModel> reaction
(
    combustionModels::psiCombustionModel::New(mesh)
);

psiReactionThermo& thermo = reaction->thermo();
thermo.validate(args.executable(), "h", "e");

volScalarField& qc = thermo.qc();
basicMultiComponentMixture& composition = thermo.composition();

//- Electron mass (default in [kg])
const dimensionedScalar eMass = constant::atomic::me;

//- Elementary charge (default in [C])
const dimensionedScalar eCharge = constant::electromagnetic::e;

//- Avogadro number (default in [1/mol])
const dimensionedScalar NA = constant::physicoChemical::NA;

//- Universal gas constant (default in [J/mol/K])
const dimensionedScalar R = constant::physicoChemical::R;

//- Boltzmann constant
const dimensionedScalar kB = constant::physicoChemical::k;

//- Pi
const scalar pi = constant::mathematical::pi;

PtrList<volScalarField>& Y = composition.Y();
const PtrList<volScalarField>& D = composition.D();

word inertSpecie(thermo.lookup("inertSpecie"));
word electronSpecie("E-");

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho()
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
const volScalarField& T = thermo.T();

#include "compressibleCreatePhi.H"

surfaceScalarField q
(
    IOobject
    (
        "q",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    linearInterpolate(U) & mesh.Sf()
);

Info << "Creating turbulence model.\n" << nl;
autoPtr<compressible::turbulenceModel> turbulence
(
    compressible::turbulenceModel::New
    (
        rho,
        U,
        phi,
        thermo
    )
);

// Set the turbulence into the reaction model
reaction->setTurbulence(turbulence());


Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
    IOobject
    (
        "dpdt",
        runTime.timeName(),
        mesh
    ),
    mesh,
    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);

Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));


multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;

forAll(Y, i)
{
    fields.add(Y[i]);
}
fields.add(thermo.he());

volScalarField dQ
(
    IOobject
    (
        "dQ",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
);

Info<< "Creating field electron number density\n" << endl;
volScalarField ne
(
    IOobject
    (
        "ne",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


Info<< "Creating field charge density\n" << endl;
volScalarField rhoq
(
    IOobject
    (
        "rhoq",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
    // (rho * qc) - (eCharge * ne)
);

Info<< "Creating field gas number density\n" << endl;
volScalarField ng("ng", p / R / T * NA);

Info<< "Creating field reduced electric field\n" << endl;
volScalarField En ("En", mag(E) / (ng));
scalarField EnTd(En.internalField() * EnToTableUnit);

Info<< "Creating field electron mobility\n" << endl;
volScalarField mue
(
    IOobject
    (
        "mue",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Creating field electron diffusivity\n" << endl;
volScalarField De
(
    IOobject
    (
        "De",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Creating field electron temperature\n" << endl;
volScalarField Te
(
    IOobject
    (
        "Te",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


/*
BE_IX::Bolos bolos;

bolos.readCollisions("./LXCat-June2013.txt");

bolos.presolve();
{
    bolos.maxwellian(2);
    forAll(rho, celli)
    {
        De[celli] = bolos.diffusivity();
        mue[celli] = bolos.mobility();
    }
    De.correctBoundaryConditions();
    mue.correctBoundaryConditions();
}
*/

Info<< "Calculating face flux field ve\n" << endl;
volVectorField Udrift
(
    IOobject
    (
        "Udrift",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    - (mue*E/ng)
);

volVectorField Uthermal
(
    IOobject
    (
        "Uthermal",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    - ((De/ng/Te)*fvc::grad(Te))
);

surfaceScalarField ve
(
    IOobject
    (
        "ve",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(Udrift+Uthermal) & mesh.Sf()
);


surfaceScalarField phi_drift
(
    IOobject
    (
        "phi_drift",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::NO_WRITE
    ),
    phi
);


surfaceScalarField phi_neutral
(
    IOobject
    (
        "phi_neutral",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::NO_WRITE
    ),
    phi
);


PtrList<surfaceScalarField> phis (composition.species().size());

forAll (composition.species(), isp)
{
    const scalar z(composition.z(isp));
    const label nCharge(z);

    if (nCharge != 0)
    {
        phis.set(isp,
            new surfaceScalarField
            (
                IOobject
                (
                    "phi." + composition.species()[isp],
                    runTime.timeName(),
                    mesh,
                    IOobject::READ_IF_PRESENT,
                    IOobject::AUTO_WRITE
                ),
                phi
            )
        );
    }
}


// plasmaWallFluxes
Info<< "Reading plasma wall flux bc control\n" << endl;

// electron wall flux
dictionary wallElectronFlux
(
    physicalProperties.subDict("wallElectronFlux")
);

word TeName(wallElectronFlux.lookup("TeName"));

wordList wallPatcheNames (wallElectronFlux.lookup("wallPatches"));
labelList wallPatcheIDs (wallPatcheNames.size(), 0);
scalarList wallReflexes (wallElectronFlux.lookup("wallReflexes"));

forAll (wallPatcheNames, pi)
{
    const word patchName = wallPatcheNames[pi];
    wallPatcheIDs[pi] 
        = mesh.boundaryMesh().findPatchID(patchName);
}

Info<< "plasma walls are \n" << wallPatcheNames << endl;

// ion wall flux

dictionary wallIonFluxes
(
    physicalProperties.subDict("wallIonFluxes")
);

const label nMaxTargets = 5;
const hashedWordList ions(wordList(wallIonFluxes.lookup("ions")));
wordList neutrals_ (ions.size()*nMaxTargets);
PtrList<hashedWordList> targetList (ions.size());
PtrList<scalarList> reflexes (ions.size());
PtrList<surfaceScalarField::GeometricBoundaryField> ionFluxBFs (ions.size());

Info<< ions.size() << " ions are \n" << ions << endl;

label nTargets = 0;

forAll (ions, iidx)
{
    const dictionary &wallIonFlux = wallIonFluxes.subDict(ions[iidx]);

    const Switch noTarget (wallIonFlux.lookupOrDefault("noTarget", false));

    const wordList targets (wallIonFlux.lookupOrDefault("neutrals", wordList()));

    if (noTarget)
    {
        targetList.set(iidx, new hashedWordList());
    }
    else if (targets.empty())
    {
        word neutralName(wallIonFlux.lookup("neutral"));
        neutrals_[nTargets] = neutralName;
        nTargets += 1;

        targetList.set(iidx, new hashedWordList(wordList(1, neutralName)));
    }
    else
    {
        wordList::subList sub (neutrals_, targets.size(), nTargets);
        forAll (targets, tidx)
        {
            sub[tidx] = targets[tidx];
        }
        nTargets += targets.size();;

        targetList.set(iidx, new hashedWordList(targets));
    }

    reflexes.set(iidx, new scalarList(wallIonFlux.lookup("wallReflexes")));

    ionFluxBFs.set(iidx,
            new surfaceScalarField::GeometricBoundaryField(phi.boundaryField()));
}

hashedWordList neutrals;

label nNeutrals = 0;

for (label nidx = 0; nidx < nTargets; nidx++)
{
    const word n(neutrals_[nidx]);

    if (!neutrals.contains(n))
    {
        neutrals.append(n);
        nNeutrals += 1;
    }
}

PtrList<surfaceScalarField::GeometricBoundaryField> neutralFluxBFs (neutrals.size());

forAll (neutrals, iidx)
{
    neutralFluxBFs.set(iidx,
            new surfaceScalarField::GeometricBoundaryField (phi.boundaryField()));
}


Info<< "Reading physicalProperties/crossSections\n" << endl;
dictionary crossSections
(
    physicalProperties.subDict("crossSections")
);

wordList csSpecies((crossSections.lookup("species")));
Info<< "Cross sections are considered for following species\n" << endl;
Info<< csSpecies << endl;

labelList csSpeciesI(csSpecies.size(), -1);

PtrList< DataEntry< scalar > > csList (csSpecies.size());

forAll (csSpecies, isp)
{
    csSpeciesI[isp] = composition.species()[csSpecies[isp]];

    csList.set(
        isp,
        DataEntry<scalar>::New(csSpecies[isp], crossSections)
    );
}

Info<< csSpeciesI << endl;

scalarField csSpeciesW (csSpecies.size(), 0.0);
forAll (csSpecies, isp)
{
    csSpeciesW[isp] = composition.W(csSpeciesI[isp]);
}

Info<< csSpeciesW << endl;


scalarField sigma_m (Te);


if (constETransport)
{
        mue.internalField() = mueNFac * mueN.value(0.0);
        mue.correctBoundaryConditions();
}