Where to start

Before any coding, define what you need the new feature to do. To detail it as possible, help yourself by building a pseudo-code in paper along with a flow-chart. Answer questions like: if you need to read probes in the mesh or directly inject numbers to a certain field, how many variables are often shared in common within work-steps, if there are repetitive tasks, and, what can help you to inspect performance while the simulation is running.

Copy to the same path the folder of the source term of openfoam that shares more similarity to your purpose. In my case, I wanted to make new sub-models of actuator disks to simulate wind turbines so, as a starting point I made a copy actuationDiskSource folder and rename this and inner files with 'calibratedActuatorDisk'. Additionally, I replaced inside each file wherever the keyword 'actuationDiskSource' was with 'calibratedActuatorDisk'. 

Note: The present content is based on OpenFOAM ESI versions, which may differ from Foundation versions in the name of some classes, variables and some constructor templates. Current syntax might raise compilation errors on that source type. 

The calibratedActuatorDisk

What is being proposed in this new source is mainly access to different methods to compute the magnitude of the thrust force, and applied it in two types of distribution. To quantify the thrust force, the variable forceMethod_ can be any of the following options: constantParameters, calafMethod, constantCt, and interpolatedFromTable. The first one is similar to actuationDiskSource:calcFroudeMethod but only depends on Ct_inf, instead of the induction factor. The next three methods requires to monitor values in order to compute the thrust force: calafMethod uses the local cellSet if the actuatorDisk, but methods constantCt, and interpolatedFromTable can read field values from an upstream point or a second cellSet. The last method requires an input file in format .csv and the execution of another user-built class called 'myLookupTable4Turbines' to interpolate the reference wind speed along with Ct_inf and Cp_inf. 

Key points about the structure:

Inheritance:

If you are not going to use the methods nor attributes that are currently implemented in the class where you copy from, then, follow my starting point instead on deriving a child from this class used as reference. An example of child-class is rotorDiskSource. As I wanted to handle different variable names and methods, I let my new fvOption class to only inheritate from cellSetOptions and writeFile classes.

Global and local variables:

Identify what variables are called within methods and external functions, and make them global to ensure their availability inside the whole class. Try to recycle variables as possible, this can be achieved by reassigning values if needed and trying to stick the algorithm to generalized or virtual functions. If a variable is used as a counter or auxiliar storing, mostly limited in the scope of a certain function, that one is local only. 

⚠️ Be cautious when using static variables, depending on the data-type and usage, it may generate conflict if you want to run more than one instance of the class (more than one source-term of the same type). In the case of conflict, your class may compile without errors but you will see unwanted performance of your source terms.  

Pack repeatitve tasks into functions and short the length of your files:

This will ease the scalability of your class and save space of lines, reducing probabilities of error when modifications are applied. You can here either use void or any template type function and try to generalize the tasks as possible. Examples of this practice in the built class: getting geometrical features of the cellSet, monitoring values of certain fields on mesh, applying source term and, writing data to outputfile, all of them are often used in any case inside the class.

Ensure paralelism using reduce():

Once you have compiled without errors and run a test case in serial mode, you may find everything working as expected, but it might not happen if you run it in parallel, some variables get crazy values and even simulation may diverge because a floating point issue. Why any of this could happen? because when our mesh is partitioned into fractions per core, some tasks may be shared within two or more processors and if we do not take care of handling how they are binding together the information, our variables may not be rightfully updated because of poor communication treatment between processors. 

 Examples of this practice in the built class: 

1. After any block of for(){}, apply reduce() in the variable you have been iterating with.  Here I used the scalar variable 'magSource_X' as a storage of the total source term magnitude that has been applied across the cellSet, to ensure we gather the TOTAL magnitude applied we use sumOp<>:

reduce(magSource_X, sumOp<scalar>());

2. As a priori, it is unknown where a given point is located inside the partitioned mesh. You can find the closest cell center to the point by comparing which processor did not return -1 for their cells, the cellID (positive integer value) will be the one assigned to 'globalCelli'. 

label localCelli = mesh_.findCell(monitorPoint_);

label globalCelli = returnReduce(localCelli, maxOp<label>());

Handle situations of error and provide information:

You can be a step forward on a floating-point error situation if you set a mathematical operation of division, add a small tolerance to the dominator in order to prevent a division by zero. 

errorRatioT = (1.0 - (Thrust_  / (magSource_X + 0.0001)));

For each situation of error in the run, manage to print in screen a message of error to help to identify the source of the issue:

1. If a point coordinate is not inside the mesh:

if (globalCelli == -1)

{

FatalErrorInFunction

<< "No owner cell found for point - Verify input coordinates: "

<< monitorPoint_ << endl;

}

2. The cellSet to work with has 0 numbers of cells:

label szMonitorCells = monitorCellSetIDs_.size();

reduce(szMonitorCells, sumOp<label>());

if (szMonitorCells == 0)

{

FatalErrorInFunction

<< "No cell is available for incoming velocity monitoring in the cellSet: "

<< monitorCellSetName_

<< exit(FatalError);

}

3. If a scalar variable has been forgotten in the dictionary:

if (Ct_star_ < VSMALL)

{

FatalErrorInFunction

<< "Scaled thrust coefficient is equal or below zero"

<< "Please check dictionary "

<< exit(FatalError);

}

4. If a string variable has been forgotten in the dictionary:

if (!csvFile_.empty())

else

{

FatalErrorInFunction

<< "-   pathToCsvTable entered not valid (check path, format and content)."

<< exit(FatalError);

}

5. If a function could not succesfully executed:

if(!lookup.findByU(Umonitor_, Ct_ref_, Cp_ref_))

{

Foam::Info << "No valid interpolation found for Umonitor_ = " << Umonitor_ << Foam::nl;

Foam::Info << "Actuatordisk: " << this->name() << " turned off." << Foam::nl;

Uref_ = 0.0;

Ct_ref_ = 0.0;

Cp_ref_ = 0.0;

}

Compile it

Go to src/fvOptions/Make/files and add the relative path to the .C file which carries the constructor. Go back to src/fvOptions and run:

wclean .

wmake

How to make your feature available in any OpenFOAM version to run in any machine here. 

You can take a look at my repo to get more details of the workflow I followed.

References:

Marc Calaf, Charles Meneveau, Johan Meyers; Large eddy simulation study of fully developed wind-turbine array boundary layers. Physics of Fluids 1 22 (1): 015110. doi: 10.1063/1.3291077

Réthoré, P.-E., van der Laan, P., Troldborg, N., Zahle, F. and Sørensen, N.N. (2014), Verification and validation of an actuator disc model. Wind Energ., 17: 919-937. doi: 10.1002/we.1607

van der Laan, M. P., Sørensen, N. N., Réthoré, P.-E., Mann, J.,  Kelly, M. C., Troldborg, N., Schepers, J. G., and Machefaux, E. (2015), An improved k- ϵ model applied to a wind turbine wake in atmospheric turbulence. Wind Energ., 18: 889–907. doi: 10.1002/we.1736.

van der Laan, M. P., Sørensen, N. N., Réthoré, P.-E., Mann, J., Kelly, M. C., and Troldborg, N. (2015) The k-ε-fP model applied to double wind turbine wakes using different actuator disk force methods. Wind Energ., 18: 2223–2240. doi: 10.1002/we.1816.

Gonzalo P. Navarro Diaz, Matias Avila and Arnau Folch (2018), An annual energy production estimation methodology for onshore wind farms over complex terrain using a RANS model with actuator discs. Journal of Physics: Conference Series, Volume 1037, Issue 7. doi: 10.1088/1742-6596/1037/7/072018