| 352 | } // setDataOnPatchHierarchy |
| 353 | |
| 354 | void |
| 355 | AdvDiffStochasticForcing::setDataOnPatch(const int data_idx, |
| 356 | Pointer<Variable<NDIM>> /*var*/, |
| 357 | Pointer<Patch<NDIM>> patch, |
| 358 | const double /*data_time*/, |
| 359 | const bool initial_time, |
| 360 | Pointer<PatchLevel<NDIM>> /*patch_level*/) |
| 361 | { |
| 362 | Pointer<CellData<NDIM, double>> divF_cc_data = patch->getPatchData(data_idx); |
| 363 | divF_cc_data->fillAll(0.0); |
| 364 | if (initial_time) return; |
| 365 | const Box<NDIM>& patch_box = patch->getBox(); |
| 366 | const Pointer<CartesianPatchGeometry<NDIM>> pgeom = patch->getPatchGeometry(); |
| 367 | const double* const dx = pgeom->getDx(); |
| 368 | double dV = 1.0; |
| 369 | for (unsigned int d = 0; d < NDIM; ++d) dV *= dx[d]; |
| 370 | const double kappa = d_adv_diff_solver->getDiffusionCoefficient(d_C_var); |
| 371 | const double dt = d_adv_diff_solver->getCurrentTimeStepSize(); |
| 372 | const double scale = d_std * std::sqrt(2.0 * kappa / (dt * dV)); |
| 373 | double C; |
| 374 | d_f_parser.DefineVar("c", &C); |
| 375 | d_f_parser.DefineVar("C", &C); |
| 376 | Pointer<CellData<NDIM, double>> C_current_cc_data = patch->getPatchData(d_C_current_cc_idx); |
| 377 | Pointer<CellData<NDIM, double>> C_half_cc_data = patch->getPatchData(d_C_half_cc_idx); |
| 378 | Pointer<CellData<NDIM, double>> C_new_cc_data = patch->getPatchData(d_C_new_cc_idx); |
| 379 | Pointer<SideData<NDIM, double>> F_sc_data = patch->getPatchData(d_F_sc_idx); |
| 380 | Pointer<CellDataFactory<NDIM, double>> C_factory = d_C_var->getPatchDataFactory(); |
| 381 | const int C_depth = C_factory->getDefaultDepth(); |
| 382 | SideData<NDIM, double> f_scale_sc_data(patch_box, C_depth, IntVector<NDIM>(0)); |
| 383 | const TimeSteppingType convective_time_stepping_type = d_adv_diff_solver->getConvectiveTimeSteppingType(d_C_var); |
| 384 | const int cycle_num = d_adv_diff_solver->getCurrentCycleNumber(); |
| 385 | for (int d = 0; d < C_depth; ++d) |
| 386 | { |
| 387 | for (int axis = 0; axis < NDIM; ++axis) |
| 388 | { |
| 389 | for (BoxIterator<NDIM> i(SideGeometry<NDIM>::toSideBox(patch_box, axis)); i; i++) |
| 390 | { |
| 391 | const hier::Index<NDIM>& ic = i(); |
| 392 | hier::Index<NDIM> ic_lower(ic); |
| 393 | ic_lower(axis) -= 1; |
| 394 | SideIndex<NDIM> is(ic, axis, SideIndex<NDIM>::Lower); |
| 395 | double f; |
| 396 | switch (convective_time_stepping_type) |
| 397 | { |
| 398 | case FORWARD_EULER: |
| 399 | { |
| 400 | C = 0.5 * ((*C_current_cc_data)(ic, d) + (*C_current_cc_data)(ic_lower, d)); |
| 401 | f = d_f_parser.Eval(); |
| 402 | break; |
| 403 | } |
| 404 | case MIDPOINT_RULE: |
| 405 | { |
| 406 | C = 0.5 * ((*C_half_cc_data)(ic, d) + (*C_half_cc_data)(ic_lower, d)); |
| 407 | f = d_f_parser.Eval(); |
| 408 | break; |
| 409 | } |
| 410 | case TRAPEZOIDAL_RULE: |
| 411 | { |
no test coverage detected