/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Website:  https://openfoam.org                  |
|   \\  /    A nd           | Version:  9                                     |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     rheoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         1;

deltaT          0.00002;

writeControl    runTime;

writeInterval   0.2;

purgeWrite      0;

writeFormat     ascii;

writePrecision  12;

writeCompression compressed;

timeFormat      general;

timePrecision   10;

graphFormat     raw;

runTimeModifiable yes;

adjustTimeStep off;

maxCo          0.1;

maxDeltaT      0.01;

functions
{

// Monitor for convergence

    probes
    {
        // Where to load it from
        functionObjectLibs ( "libsampling.so" );

        type            probes;

        // Name of the directory for probe data
        name            probes;

        // Write at same frequency as fields
        writeControl   timeStep;
        writeInterval  50;

        // Fields to be probed
        fields
        (
            U
	    p
	    tau
        );

        probeLocations
        (
             (0.001 0.000731571 0.000731571)
	     (0.00592628 0.000731571 0.000731571)
        );
    }
	
    outputCd
    {
	 functionObjectLibs ("libutilityFunctionObjects.so");
	 type coded;
	 name calCd;
	 writeControl timeStep;
         writeInterval      10;
         enabled             yes;

	 codeWrite
	 #{
 
          // Lookup/create variable 
	     
	   const volVectorField& U = mesh().lookupObject<volVectorField>("U");
           const volSymmTensorField& tau = mesh().lookupObject<volSymmTensorField>("tau");
           const volScalarField& p = mesh().lookupObject<volScalarField>("p");
           const dictionary& constDict = mesh().lookupObject<IOdictionary>("constitutiveProperties");
           dimensionedScalar rho_(constDict.subDict("parameters").lookup("rho"));
           dimensionedScalar etaS_(constDict.subDict("parameters").lookup("etaS"));
           dimensionedScalar etaP_(constDict.subDict("parameters").lookup("etaP"));
  
           label sph = mesh().boundaryMesh().findPatchID("Sphere3-5");
           scalarList list;
    
          // Compute cd
 
           volTensorField L(fvc::grad(U));

	   volSymmTensorField F(tau + symm( L + L.T() ) * etaS_ - p * symmTensor::I * rho_);

           vector Fpatch = gSum( ( -mesh().boundaryMesh()[sph].faceAreas() ) & F.boundaryField()[sph] );
  
           list.append(mesh().time().value());  // Time (col 0)  
           list.append(Fpatch.x());             // F   (col 1)  

          // Write data

           string comsh;           
           string filename("force.txt");
	   std::stringstream doub2str; doub2str.precision(12);

           comsh = "./writeData " + filename;
           forAll(list, id)
            {
              doub2str.str(std::string());
              doub2str << list[id]; 
              comsh += " " + doub2str.str();
            }
           
           if (Pstream::master())
            {
	      system(comsh);
            }

	 #};
    }

    QField
    {
        type coded;
        libs            ("libutilityFunctionObjects.so");
        writeControl    writeTime;
        executeControl  writeTime;

        codeExecute
        #{
            const volVectorField& U = mesh().lookupObject<volVectorField>("U");

            volTensorField gradU(fvc::grad(U));
            volSymmTensorField D(symm(gradU));
            volTensorField W(0.5*(gradU - gradU.T()));

            volScalarField DD(magSqr(D));
            volScalarField WW(magSqr(W));

            volScalarField Q
            (
                IOobject
                (
                    "Q",
                    mesh().time().timeName(),
                    mesh(),
                    IOobject::NO_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                (DD - WW)/(DD + WW + dimensionedScalar("tiny", DD.dimensions(), SMALL))
            );
            Q.write();
        #};
    }

    #includeFunc residuals
}
// ************************************************************************* //
