Actual source code: dtprob.c

      1: #include <petscdt.h>

      3: #include <petscvec.h>
      4: #include <petscdraw.h>
      5: #include <petsc/private/petscimpl.h>

      7: const char *const DTProbDensityTypes[] = {"constant", "gaussian", "maxwell_boltzmann", "DTProb Density", "DTPROB_DENSITY_", NULL};

      9: /*@
     10:   PetscPDFMaxwellBoltzmann1D - PDF for the Maxwell-Boltzmann distribution in 1D

     12:   Not Collective

     14:   Input Parameters:
     15: + x     - Speed in $[0, \infty]$
     16: - dummy - Unused

     18:   Output Parameter:
     19: . p - The probability density at `x`

     21:   Level: beginner

     23: .seealso: `PetscCDFMaxwellBoltzmann1D()`
     24: @*/
     25: PetscErrorCode PetscPDFMaxwellBoltzmann1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
     26: {
     27:   p[0] = PetscSqrtReal(2. / PETSC_PI) * PetscExpReal(-0.5 * PetscSqr(x[0]));
     28:   return PETSC_SUCCESS;
     29: }

     31: /*@
     32:   PetscCDFMaxwellBoltzmann1D - CDF for the Maxwell-Boltzmann distribution in 1D

     34:   Not Collective

     36:   Input Parameters:
     37: + x     - Speed in $[0, \infty]$
     38: - dummy - Unused

     40:   Output Parameter:
     41: . p - The probability density at `x`

     43:   Level: beginner

     45: .seealso: `PetscPDFMaxwellBoltzmann1D()`
     46: @*/
     47: PetscErrorCode PetscCDFMaxwellBoltzmann1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
     48: {
     49:   p[0] = PetscErfReal(x[0] / PETSC_SQRT2);
     50:   return PETSC_SUCCESS;
     51: }

     53: /*@
     54:   PetscPDFMaxwellBoltzmann2D - PDF for the Maxwell-Boltzmann distribution in 2D

     56:   Not Collective

     58:   Input Parameters:
     59: + x     - Speed in $[0, \infty]$
     60: - dummy - Unused

     62:   Output Parameter:
     63: . p - The probability density at `x`

     65:   Level: beginner

     67: .seealso: `PetscCDFMaxwellBoltzmann2D()`
     68: @*/
     69: PetscErrorCode PetscPDFMaxwellBoltzmann2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
     70: {
     71:   p[0] = x[0] * PetscExpReal(-0.5 * PetscSqr(x[0]));
     72:   return PETSC_SUCCESS;
     73: }

     75: /*@
     76:   PetscCDFMaxwellBoltzmann2D - CDF for the Maxwell-Boltzmann distribution in 2D

     78:   Not Collective

     80:   Input Parameters:
     81: + x     - Speed in $[0, \infty]$
     82: - dummy - Unused

     84:   Output Parameter:
     85: . p - The probability density at `x`

     87:   Level: beginner

     89: .seealso: `PetscPDFMaxwellBoltzmann2D()`
     90: @*/
     91: PetscErrorCode PetscCDFMaxwellBoltzmann2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
     92: {
     93:   p[0] = 1. - PetscExpReal(-0.5 * PetscSqr(x[0]));
     94:   return PETSC_SUCCESS;
     95: }

     97: /*@
     98:   PetscPDFMaxwellBoltzmann3D - PDF for the Maxwell-Boltzmann distribution in 3D

    100:   Not Collective

    102:   Input Parameters:
    103: + x     - Speed in $[0, \infty]$
    104: - dummy - Unused

    106:   Output Parameter:
    107: . p - The probability density at `x`

    109:   Level: beginner

    111: .seealso: `PetscCDFMaxwellBoltzmann3D()`
    112: @*/
    113: PetscErrorCode PetscPDFMaxwellBoltzmann3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    114: {
    115:   p[0] = PetscSqrtReal(2. / PETSC_PI) * PetscSqr(x[0]) * PetscExpReal(-0.5 * PetscSqr(x[0]));
    116:   return PETSC_SUCCESS;
    117: }

    119: /*@
    120:   PetscCDFMaxwellBoltzmann3D - CDF for the Maxwell-Boltzmann distribution in 3D

    122:   Not Collective

    124:   Input Parameters:
    125: + x     - Speed in $[0, \infty]$
    126: - dummy - Unused

    128:   Output Parameter:
    129: . p - The probability density at `x`

    131:   Level: beginner

    133: .seealso: `PetscPDFMaxwellBoltzmann3D()`
    134: @*/
    135: PetscErrorCode PetscCDFMaxwellBoltzmann3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    136: {
    137:   p[0] = PetscErfReal(x[0] / PETSC_SQRT2) - PetscSqrtReal(2. / PETSC_PI) * x[0] * PetscExpReal(-0.5 * PetscSqr(x[0]));
    138:   return PETSC_SUCCESS;
    139: }

    141: /*@
    142:   PetscPDFGaussian1D - PDF for the Gaussian distribution in 1D

    144:   Not Collective

    146:   Input Parameters:
    147: + x     - Coordinate in $[-\infty, \infty]$
    148: - scale - Scaling value

    150:   Output Parameter:
    151: . p - The probability density at `x`

    153:   Level: beginner

    155: .seealso: `PetscPDFMaxwellBoltzmann3D()`
    156: @*/
    157: PetscErrorCode PetscPDFGaussian1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
    158: {
    159:   const PetscReal sigma = scale ? scale[0] : 1.;
    160:   p[0]                  = PetscSqrtReal(1. / (2. * PETSC_PI)) * PetscExpReal(-0.5 * PetscSqr(x[0] / sigma)) / sigma;
    161:   return PETSC_SUCCESS;
    162: }

    164: PetscErrorCode PetscCDFGaussian1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
    165: {
    166:   const PetscReal sigma = scale ? scale[0] : 1.;
    167:   p[0]                  = 0.5 * (1. + PetscErfReal(x[0] / PETSC_SQRT2 / sigma));
    168:   return PETSC_SUCCESS;
    169: }

    171: /*@
    172:   PetscPDFSampleGaussian1D - Sample uniformly from a Gaussian distribution in 1D

    174:   Not Collective

    176:   Input Parameters:
    177: + p     - A uniform variable on $[0, 1]$
    178: - dummy - ignored

    180:   Output Parameter:
    181: . x - Coordinate in $[-\infty, \infty]$

    183:   Level: beginner

    185:   Note:
    186:   See <http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function> and
    187:   <https://stackoverflow.com/questions/27229371/inverse-error-function-in-c>

    189: .seealso: `PetscPDFGaussian2D()`
    190: @*/
    191: PetscErrorCode PetscPDFSampleGaussian1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
    192: {
    193:   const PetscReal q       = 2 * p[0] - 1.;
    194:   const PetscInt  maxIter = 100;
    195:   PetscReal       ck[100], r = 0.;
    196:   PetscInt        k, m;

    198:   PetscFunctionBeginHot;
    199:   /* Transform input to [-1, 1] since the code below computes the inverse error function */
    200:   for (k = 0; k < maxIter; ++k) ck[k] = 0.;
    201:   ck[0] = 1;
    202:   r     = ck[0] * (PetscSqrtReal(PETSC_PI) / 2.) * q;
    203:   for (k = 1; k < maxIter; ++k) {
    204:     const PetscReal temp = 2. * k + 1.;

    206:     for (m = 0; m <= k - 1; ++m) {
    207:       PetscReal denom = (m + 1.) * (2. * m + 1.);

    209:       ck[k] += (ck[m] * ck[k - 1 - m]) / denom;
    210:     }
    211:     r += (ck[k] / temp) * PetscPowReal((PetscSqrtReal(PETSC_PI) / 2.) * q, 2. * k + 1.);
    212:   }
    213:   /* Scale erfinv() by \sqrt{\pi/2} */
    214:   x[0] = PetscSqrtReal(PETSC_PI * 0.5) * r;
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: /*@
    219:   PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D

    221:   Not Collective

    223:   Input Parameters:
    224: + x     - Coordinate in $[-\infty, \infty]^2$
    225: - dummy - ignored

    227:   Output Parameter:
    228: . p - The probability density at `x`

    230:   Level: beginner

    232: .seealso: `PetscPDFSampleGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
    233: @*/
    234: PetscErrorCode PetscPDFGaussian2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    235: {
    236:   p[0] = (1. / PETSC_PI) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1])));
    237:   return PETSC_SUCCESS;
    238: }

    240: /*@
    241:   PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D
    242:   <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>

    244:   Not Collective

    246:   Input Parameters:
    247: + p     - A uniform variable on $[0, 1]^2$
    248: - dummy - ignored

    250:   Output Parameter:
    251: . x - Coordinate in $[-\infty, \infty]^2 $

    253:   Level: beginner

    255: .seealso: `PetscPDFGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
    256: @*/
    257: PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
    258: {
    259:   const PetscReal mag = PetscSqrtReal(-2.0 * PetscLogReal(p[0]));
    260:   x[0]                = mag * PetscCosReal(2.0 * PETSC_PI * p[1]);
    261:   x[1]                = mag * PetscSinReal(2.0 * PETSC_PI * p[1]);
    262:   return PETSC_SUCCESS;
    263: }

    265: /*@
    266:   PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D

    268:   Not Collective

    270:   Input Parameters:
    271: + x     - Coordinate in $[-\infty, \infty]^3$
    272: - dummy - ignored

    274:   Output Parameter:
    275: . p - The probability density at `x`

    277:   Level: beginner

    279: .seealso: `PetscPDFSampleGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
    280: @*/
    281: PetscErrorCode PetscPDFGaussian3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    282: {
    283:   p[0] = (1. / PETSC_PI * PetscSqrtReal(PETSC_PI)) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]) + PetscSqr(x[2])));
    284:   return PETSC_SUCCESS;
    285: }

    287: /*@
    288:   PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D
    289:   <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>

    291:   Not Collective

    293:   Input Parameters:
    294: + p     - A uniform variable on $[0, 1]^3$
    295: - dummy - ignored

    297:   Output Parameter:
    298: . x - Coordinate in $[-\infty, \infty]^3$

    300:   Level: beginner

    302: .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
    303: @*/
    304: PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
    305: {
    306:   PetscCall(PetscPDFSampleGaussian1D(p, dummy, x));
    307:   PetscCall(PetscPDFSampleGaussian2D(&p[1], dummy, &x[1]));
    308:   return PETSC_SUCCESS;
    309: }

    311: /*@
    312:   PetscPDFConstant1D - PDF for the uniform distribution in 1D

    314:   Not Collective

    316:   Input Parameters:
    317: + x     - Coordinate in $[-1, 1]$
    318: - dummy - Unused

    320:   Output Parameter:
    321: . p - The probability density at `x`

    323:   Level: beginner

    325: .seealso: `PetscCDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscPDFConstant2D()`, `PetscPDFConstant3D()`
    326: @*/
    327: PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    328: {
    329:   p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.;
    330:   return PETSC_SUCCESS;
    331: }

    333: /*@
    334:   PetscCDFConstant1D - CDF for the uniform distribution in 1D

    336:   Not Collective

    338:   Input Parameters:
    339: + x     - Coordinate in $[-1, 1]$
    340: - dummy - Unused

    342:   Output Parameter:
    343: . p - The cumulative probability at `x`

    345:   Level: beginner

    347: .seealso: `PetscPDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscCDFConstant2D()`, `PetscCDFConstant3D()`
    348: @*/
    349: PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    350: {
    351:   p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.));
    352:   return PETSC_SUCCESS;
    353: }

    355: /*@
    356:   PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D

    358:   Not Collective

    360:   Input Parameters:
    361: + p     - A uniform variable on $[0, 1]$
    362: - dummy - Unused

    364:   Output Parameter:
    365: . x - Coordinate in $[-1, 1]$

    367:   Level: beginner

    369: .seealso: `PetscPDFConstant1D()`, `PetscCDFConstant1D()`, `PetscPDFSampleConstant2D()`, `PetscPDFSampleConstant3D()`
    370: @*/
    371: PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
    372: {
    373:   x[0] = 2. * p[0] - 1.;
    374:   return PETSC_SUCCESS;
    375: }

    377: /*@
    378:   PetscPDFConstant2D - PDF for the uniform distribution in 2D

    380:   Not Collective

    382:   Input Parameters:
    383: + x     - Coordinate in $[-1, 1]^2$
    384: - dummy - Unused

    386:   Output Parameter:
    387: . p - The probability density at `x`

    389:   Level: beginner

    391: .seealso: `PetscCDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscPDFConstant1D()`, `PetscPDFConstant3D()`
    392: @*/
    393: PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    394: {
    395:   p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.;
    396:   return PETSC_SUCCESS;
    397: }

    399: /*@
    400:   PetscCDFConstant2D - CDF for the uniform distribution in 2D

    402:   Not Collective

    404:   Input Parameters:
    405: + x     - Coordinate in $[-1, 1]^2$
    406: - dummy - Unused

    408:   Output Parameter:
    409: . p - The cumulative probability at `x`

    411:   Level: beginner

    413: .seealso: `PetscPDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscCDFConstant1D()`, `PetscCDFConstant3D()`
    414: @*/
    415: PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    416: {
    417:   p[0] = x[0] <= -1. || x[1] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)) * (x[1] > 1. ? 1. : 0.5 * (x[1] + 1.));
    418:   return PETSC_SUCCESS;
    419: }

    421: /*@
    422:   PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on $[-1, 1]^2$ in 2D

    424:   Not Collective

    426:   Input Parameters:
    427: + p     - Two uniform variables on $[0, 1]$
    428: - dummy - Unused

    430:   Output Parameter:
    431: . x - Coordinate in $[-1, 1]^2$

    433:   Level: beginner

    435: .seealso: `PetscPDFConstant2D()`, `PetscCDFConstant2D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant3D()`
    436: @*/
    437: PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
    438: {
    439:   x[0] = 2. * p[0] - 1.;
    440:   x[1] = 2. * p[1] - 1.;
    441:   return PETSC_SUCCESS;
    442: }

    444: /*@
    445:   PetscPDFConstant3D - PDF for the uniform distribution in 3D

    447:   Not Collective

    449:   Input Parameters:
    450: + x     - Coordinate in $[-1, 1]^3$
    451: - dummy - Unused

    453:   Output Parameter:
    454: . p - The probability density at `x`

    456:   Level: beginner

    458: .seealso: `PetscCDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
    459: @*/
    460: PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    461: {
    462:   p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.;
    463:   return PETSC_SUCCESS;
    464: }

    466: /*@
    467:   PetscCDFConstant3D - CDF for the uniform distribution in 3D

    469:   Not Collective

    471:   Input Parameters:
    472: + x     - Coordinate in $[-1, 1]^3$
    473: - dummy - Unused

    475:   Output Parameter:
    476: . p - The cumulative probability at `x`

    478:   Level: beginner

    480: .seealso: `PetscPDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscCDFConstant1D()`, `PetscCDFConstant2D()`
    481: @*/
    482: PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
    483: {
    484:   p[0] = x[0] <= -1. || x[1] <= -1. || x[2] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)) * (x[1] > 1. ? 1. : 0.5 * (x[1] + 1.)) * (x[2] > 1. ? 1. : 0.5 * (x[2] + 1.));
    485:   return PETSC_SUCCESS;
    486: }

    488: /*@
    489:   PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on $[-1, 1]^3$ in 3D

    491:   Not Collective

    493:   Input Parameters:
    494: + p     - Three uniform variables on $[0, 1]$
    495: - dummy - Unused

    497:   Output Parameter:
    498: . x - Coordinate in $[-1, 1]^3$

    500:   Level: beginner

    502: .seealso: `PetscPDFConstant3D()`, `PetscCDFConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
    503: @*/
    504: PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
    505: {
    506:   x[0] = 2. * p[0] - 1.;
    507:   x[1] = 2. * p[1] - 1.;
    508:   x[2] = 2. * p[2] - 1.;
    509:   return PETSC_SUCCESS;
    510: }

    512: /*@C
    513:   PetscProbCreateFromOptions - Return the probability distribution specified by the arguments and options

    515:   Not Collective

    517:   Input Parameters:
    518: + dim    - The dimension of sample points
    519: . prefix - The options prefix, or `NULL`
    520: - name   - The options database name for the probability distribution type

    522:   Output Parameters:
    523: + pdf     - The PDF of this type, or `NULL`
    524: . cdf     - The CDF of this type, or `NULL`
    525: - sampler - The PDF sampler of this type, or `NULL`

    527:   Level: intermediate

    529: .seealso: `PetscProbFunc`, `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()`
    530: @*/
    531: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler)
    532: {
    533:   DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN;

    535:   PetscFunctionBegin;
    536:   PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
    537:   PetscCall(PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL));
    538:   PetscOptionsEnd();

    540:   if (pdf) {
    541:     PetscAssertPointer(pdf, 4);
    542:     *pdf = NULL;
    543:   }
    544:   if (cdf) {
    545:     PetscAssertPointer(cdf, 5);
    546:     *cdf = NULL;
    547:   }
    548:   if (sampler) {
    549:     PetscAssertPointer(sampler, 6);
    550:     *sampler = NULL;
    551:   }
    552:   switch (den) {
    553:   case DTPROB_DENSITY_CONSTANT:
    554:     switch (dim) {
    555:     case 1:
    556:       if (pdf) *pdf = PetscPDFConstant1D;
    557:       if (cdf) *cdf = PetscCDFConstant1D;
    558:       if (sampler) *sampler = PetscPDFSampleConstant1D;
    559:       break;
    560:     case 2:
    561:       if (pdf) *pdf = PetscPDFConstant2D;
    562:       if (cdf) *cdf = PetscCDFConstant2D;
    563:       if (sampler) *sampler = PetscPDFSampleConstant2D;
    564:       break;
    565:     case 3:
    566:       if (pdf) *pdf = PetscPDFConstant3D;
    567:       if (cdf) *cdf = PetscCDFConstant3D;
    568:       if (sampler) *sampler = PetscPDFSampleConstant3D;
    569:       break;
    570:     default:
    571:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
    572:     }
    573:     break;
    574:   case DTPROB_DENSITY_GAUSSIAN:
    575:     switch (dim) {
    576:     case 1:
    577:       if (pdf) *pdf = PetscPDFGaussian1D;
    578:       if (cdf) *cdf = PetscCDFGaussian1D;
    579:       if (sampler) *sampler = PetscPDFSampleGaussian1D;
    580:       break;
    581:     case 2:
    582:       if (pdf) *pdf = PetscPDFGaussian2D;
    583:       if (sampler) *sampler = PetscPDFSampleGaussian2D;
    584:       break;
    585:     case 3:
    586:       if (pdf) *pdf = PetscPDFGaussian3D;
    587:       if (sampler) *sampler = PetscPDFSampleGaussian3D;
    588:       break;
    589:     default:
    590:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
    591:     }
    592:     break;
    593:   case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
    594:     switch (dim) {
    595:     case 1:
    596:       if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
    597:       if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
    598:       break;
    599:     case 2:
    600:       if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
    601:       if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
    602:       break;
    603:     case 3:
    604:       if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
    605:       if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
    606:       break;
    607:     default:
    608:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
    609:     }
    610:     break;
    611:   default:
    612:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
    613:   }
    614:   PetscFunctionReturn(PETSC_SUCCESS);
    615: }

    617: #ifdef PETSC_HAVE_KS
    618: EXTERN_C_BEGIN
    619:   #include <KolmogorovSmirnovDist.h>
    620: EXTERN_C_END
    621: #endif

    623: /*@C
    624:   PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.

    626:   Collective

    628:   Input Parameters:
    629: + v   - The data vector, blocksize is the sample dimension
    630: - cdf - The analytic CDF

    632:   Output Parameter:
    633: . alpha - The KS statistic

    635:   Level: advanced

    637:   Notes:
    638:   The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is

    640:   $$
    641:   D_n = \sup_x \left| F_n(x) - F(x) \right|
    642:   $$

    644:   where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by

    646:   $$
    647:   F_n =  # of samples <= x / n
    648:   $$

    650:   The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
    651:   cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
    652:   the largest absolute difference between the two distribution functions across all $x$ values.

    654: .seealso: `PetscProbFunc`
    655: @*/
    656: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha)
    657: {
    658: #if !defined(PETSC_HAVE_KS)
    659:   SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
    660: #else
    661:   PetscViewer        viewer = NULL;
    662:   PetscViewerFormat  format;
    663:   PetscDraw          draw;
    664:   PetscDrawLG        lg;
    665:   PetscDrawAxis      axis;
    666:   const PetscScalar *a;
    667:   PetscReal         *speed, Dn = PETSC_MIN_REAL;
    668:   PetscBool          isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg;
    669:   PetscInt           n, p, dim, d;
    670:   PetscMPIInt        size;
    671:   const char        *names[2] = {"Analytic", "Empirical"};
    672:   char               title[PETSC_MAX_PATH_LEN];
    673:   PetscOptions       options;
    674:   const char        *prefix;
    675:   MPI_Comm           comm;

    677:   PetscFunctionBegin;
    678:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
    679:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, &prefix));
    680:   PetscCall(PetscObjectGetOptions((PetscObject)v, &options));
    681:   PetscCall(PetscOptionsCreateViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg));
    682:   if (flg) {
    683:     PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    684:     PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    685:   }
    686:   if (isascii) {
    687:     PetscCall(PetscViewerPushFormat(viewer, format));
    688:   } else if (isdraw) {
    689:     PetscCall(PetscViewerPushFormat(viewer, format));
    690:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    691:     PetscCall(PetscDrawLGCreate(draw, 2, &lg));
    692:     PetscCall(PetscDrawLGSetLegend(lg, names));
    693:   }

    695:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
    696:   PetscCall(VecGetLocalSize(v, &n));
    697:   PetscCall(VecGetBlockSize(v, &dim));
    698:   n /= dim;
    699:   /* TODO Parallel algorithm is harder */
    700:   if (size == 1) {
    701:     PetscCall(PetscMalloc1(n, &speed));
    702:     PetscCall(VecGetArrayRead(v, &a));
    703:     for (p = 0; p < n; ++p) {
    704:       PetscReal mag = 0.;

    706:       for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
    707:       speed[p] = PetscSqrtReal(mag);
    708:     }
    709:     PetscCall(PetscSortReal(n, speed));
    710:     PetscCall(VecRestoreArrayRead(v, &a));
    711:     for (p = 0; p < n; ++p) {
    712:       const PetscReal x = speed[p], Fn = ((PetscReal)p) / n;
    713:       PetscReal       F, vals[2];

    715:       PetscCall(cdf(&x, NULL, &F));
    716:       Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
    717:       if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
    718:       if (isdraw) {
    719:         vals[0] = F;
    720:         vals[1] = Fn;
    721:         PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
    722:       }
    723:     }
    724:     if (speed[n - 1] < 6.) {
    725:       const PetscReal k = (PetscInt)(6. - speed[n - 1]) + 1, dx = (6. - speed[n - 1]) / k;
    726:       for (p = 0; p <= k; ++p) {
    727:         const PetscReal x = speed[n - 1] + p * dx, Fn = 1.0;
    728:         PetscReal       F, vals[2];

    730:         PetscCall(cdf(&x, NULL, &F));
    731:         Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
    732:         if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
    733:         if (isdraw) {
    734:           vals[0] = F;
    735:           vals[1] = Fn;
    736:           PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
    737:         }
    738:       }
    739:     }
    740:     PetscCall(PetscFree(speed));
    741:   }
    742:   if (isdraw) {
    743:     PetscCall(PetscDrawLGGetAxis(lg, &axis));
    744:     PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn));
    745:     PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)"));
    746:     PetscCall(PetscDrawLGDraw(lg));
    747:     PetscCall(PetscDrawLGDestroy(&lg));
    748:   }
    749:   if (viewer) {
    750:     PetscCall(PetscViewerFlush(viewer));
    751:     PetscCall(PetscViewerPopFormat(viewer));
    752:     PetscCall(PetscViewerDestroy(&viewer));
    753:   }
    754:   *alpha = KSfbar((int)n, (double)Dn);
    755:   PetscFunctionReturn(PETSC_SUCCESS);
    756: #endif
    757: }