| 623 | |
| 624 | |
| 625 | void ca_dercircvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, |
| 626 | const FArrayBox& datfab, const Geometry& geom, |
| 627 | Real /*time*/, const int* /*bcrec*/, int /*level*/) |
| 628 | { |
| 629 | |
| 630 | // our input dat is rho, UMX, UMY, UMZ |
| 631 | |
| 632 | auto const dat = datfab.array(); |
| 633 | auto const der = derfab.array(); |
| 634 | |
| 635 | auto dx = geom.CellSizeArray(); |
| 636 | auto problo = geom.ProbLoArray(); |
| 637 | auto geomdata = geom.data(); |
| 638 | |
| 639 | amrex::ParallelFor(bx, |
| 640 | [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept |
| 641 | { |
| 642 | |
| 643 | GpuArray<Real, 3> loc = {0.0}; |
| 644 | loc[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0]; |
| 645 | #if AMREX_SPACEDIM >= 2 |
| 646 | loc[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1]; |
| 647 | #endif |
| 648 | #if AMREX_SPACEDIM == 3 |
| 649 | loc[2] = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2]; |
| 650 | #endif |
| 651 | |
| 652 | if (domain_is_plane_parallel) { |
| 653 | #if AMREX_SPACEDIM == 2 |
| 654 | // the circumferential velocity is just the out-of-plane velocity |
| 655 | der(i,j,k,0) = dat(i,j,k,3)/dat(i,j,k,0); |
| 656 | #elif AMREX_SPACEDIM == 3 |
| 657 | // the velocity in the x-y plane decomposed into r, phi unit vectors is: |
| 658 | // v_cyl = ( u cos phi + v sin phi) e_r + |
| 659 | // (-u sin phi + v cos phi) e_phi |
| 660 | // where e_r and e_phi are the cylindrical unit vectors |
| 661 | |
| 662 | // we need the distance in the x-y plane from the origin |
| 663 | Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); |
| 664 | der(i,j,k,0) = (-dat(i,j,k,1)*loc[1] + dat(i,j,k,2)*loc[0]) / (dat(i,j,k,0)*r); |
| 665 | #endif |
| 666 | } else { |
| 667 | Real r = distance(geomdata, loc); |
| 668 | |
| 669 | // we really mean just the velocity component that is |
| 670 | // perpendicular to radial, and in general 3-d (e.g. a |
| 671 | // sphere), the sign doesn't make sense, so we compute this |
| 672 | // such that v_r^2 + v_c^2 = v^2 |
| 673 | Real vtot2 = (dat(i,j,k,1)*dat(i,j,k,1) + |
| 674 | dat(i,j,k,2)*dat(i,j,k,2) + |
| 675 | dat(i,j,k,3)*dat(i,j,k,3))/(dat(i,j,k,0)*dat(i,j,k,0)); |
| 676 | |
| 677 | Real vr = (dat(i,j,k,1)*loc[0] + |
| 678 | dat(i,j,k,2)*loc[1] + |
| 679 | dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r ); |
| 680 | |
| 681 | der(i,j,k,0) = std::sqrt(amrex::max(vtot2 - vr*vr, 0.0_rt)); |
| 682 | } |
nothing calls this directly
no test coverage detected