Parameter study in 3D parameter space

Problem description

This tutorial is an extension to ParameterStudy2DDemo, and therefore has similar code structure. In ParameterStudy1DDemo and ParameterStudy2DDemo it has been shown how to do 1D and 2D parameter studies. In those cases only the particle diameter has been changed, based on the simulation that is run. In this tutorial it is shown how to perform a 3D parameter study on contact model properties, more specifically on studying the effect of stiffness, dissipation and density. The parameter space now looks as follows:

Schematic overview of a 3D parameter study. Definitions of i, j ,k , dim1, dim2 and dim3 will be given later.

Headers

The following headers are included:

#include "Mercury3D.h"
#include "Math/Helpers.h"
#include "Species/LinearViscoelasticFrictionSpecies.h"
#include <iostream>
#include <vector>

These are the headers needed for this tutorial.

Before the main function

The main class inherits from Mercury3D.

ParameterStudy3DDemo : public Mercury3D 

All steps needed in order to create this class are explained below.
Step 1: (Private) Member variables.
Similar to the extension of 1D to 2D parameter studies a new integer is introduced to store the size of the third dimension of the parameter study. Furthermore the vector studyNum is again extended with 1, as it needs to store the information of the third dimension. studyNum is now build up as follows: The first entry is (still) storing an integer value giving information of the completeness of the parameter study. The second, third and fourth entry store the current run in the first, second and third dimension respectively.

    int dim1_ = 0;
    int dim2_ = 0;
    int dim3_ = 0;
    std::vector<int> studyNum;

Step 2: Set- and Get-functions for private member variables.
The set- and get-functions are extended to incorporate the third dimension.

    void setStudyDimensions(int dim1, int dim2, int dim3)
    {
        dim1_ = dim1;
        dim2_ = dim2;
        dim3_ = dim3;
    }
 
    std::vector<int> getCurrentStudyNum()
    {
        return DPMBase::get3DParametersFromRunNumber(dim1_, dim2_, dim3_);
    }

Step 3: Create the particle species.
This tutorial was build around the wish to perform a 3D parameter study over the contact model properties. Hence, compared to ParameterStudy1DDemo and ParameterStudy2DDemo, setting the particle species is now important. The code however is still very similar, as createSpecies() had already been incorporated in the class. For this reason we can simply use the same functions as before to get the simulation run numbers. In this case it has been chosen to vary the stiffness, dissipation and density of the particles. However, if desired, all other particle properties could be changed by the run numbers.

    void createSpecies()
    {
        //Lets say we want to update the particle stiffness, dissipation and density based on studyNum.
        double k_ref = 1.0e5;
        double diss_ref = 0.63;
        double rho_ref = 2000.0;
 
        // Arbitrary change to the values based on the three study numbers stored in studyNum
        double k_run = k_ref * studyNum[1];
        double diss_run = diss_ref + 0.05*studyNum[2];
        double rho_run = rho_ref + (studyNum[3] - 1)*500;
 
        // Setup the contact model for the species
        LinearViscoelasticFrictionSpecies species;
        species.setDensity(rho_run);
        //specify contact properties
        //normal forces
        species.setStiffness(k_run);
        species.setDissipation(diss_run);
        //tangential (sliding) forces
        species.setSlidingFrictionCoefficient(0.5);
        species.setSlidingStiffness(1.2e4);
        species.setSlidingDissipation(0.16);
        //tangential (rolling) torques
        species.setRollingFrictionCoefficient(0.2);
        species.setRollingStiffness(1.2e4);
        species.setRollingDissipation(6.3e-2);
        //normal (torsion/spin) torques
        species.setTorsionFrictionCoefficient(0.1);
        species.setTorsionStiffness(1.2e4);
        species.setSlidingDissipation(6.3e-2);
        speciesHandler.copyAndAddObject(species);
    }

Step 4: Setting up the initial conditions.
The only difference with respect to ParameterStudy1DDemo and ParameterStudy2DDemo is that now the particles do not vary in size anymore.

    void setupInitialConditions() override
    {
        studyNum = getCurrentStudyNum();
        createSpecies();
 
        SphericalParticle P;
        P.setSpecies(speciesHandler.getObject(0));
        // Create 2 particles
        particleHandler.copyAndAddObject(P);
        particleHandler.copyAndAddObject(P);
        //Set up the stuff that is the same for all runs
        particleHandler.getObject(0)->setPosition(Vec3D(0.006,0.005,0.0));
        particleHandler.getObject(1)->setPosition(Vec3D(0.004,0.005,0.0));
        particleHandler.getObject(0)->setVelocity(Vec3D(-0.1,0.0,0.0));
        particleHandler.getObject(1)->setVelocity(Vec3D( 0.1,0.0,0.0));
        // We have two different size particles for both left and right particles
        // The value gets updated based on what studyNum is, these come from the automatic counter
        particleHandler.getObject(0)->setRadius(0.0001);
        particleHandler.getObject(1)->setRadius(0.0001);
    }

Step 5: Setting up the parameter study.
Nothing changed with respect to ParameterStudy2DDemo, both actionsBeforeTimeLoop() and actionsAfterSolve() are the same.

Step 6: Creating the main function.
The only difference with ParameterStudy3DDemo is that now all three study dimensions have to be set.

problem.setStudyDimensions(2,2,3) 

This line results in 12 simulations being run, in which all combinations of the contact model parameters are solved.

int main(int argc UNUSED, char *argv[] UNUSED)
{
    //Problem constructor, using class definition above
    ParameterStudy3DDemo problem;
    // Setup the study dimensions that you want to perform
    problem.setStudyDimensions(2,2,3);
    // Setup the automatic numbering system (create/read COUNTER_DO_NOT_DEL file)
    problem.autoNumber();
    // Setup output file name
    problem.setName("ParameterStudy3DDemo");
    // General simulation settings
    problem.setTimeStep(1e-5);
    problem.setTimeMax(0.1);
    problem.setSaveCount(5000);
    problem.setXMin(0.0);
    problem.setXMax(1.0);
    problem.setYMin(0.0);
    problem.setYMax(1.0);
    problem.setZMin(0.0);
    problem.setZMax(1.0);
 
    // Solve the problem
    problem.solve();
    // Give a successful exit code to the terminal
    return 0;
}

If you would like to look at the 1D or 2D parameter study demos:
ParameterStudy1DDemo
ParameterStudy2DDemo
Or return to Overview of advanced tutorials