  We model the structure and evolution of black hole accretion disks, and their
neighboring regions, using numerical simulations. The numerics is governed by
the equations of general relativistic magneto-hydrodynamics (GRMHD). In
particular, such disks and outflows can be found at the base of very energetic
ultra-relativistic jets produced by cosmic explosions, so called gamma-ray
bursts (GRBs). Another, more persistent type of the jet phenomena, are blazars,
emitted from the centers of galaxies. Long-lasting, detailed computations are
essential to properly determine the physics of these explosions, and confront
the theoretical models with any potential observables. From the point of view
of numerical methods and computational techniques, three ingredients need to be
considered. First, the numerical scheme must work in a conservative manner,
which is achieved by solving a set of non-linear equations at each time-step,
to advance the conserved quantities from one time step to the next. Second, the
efficiency of computations intrinsically depends on the code parallelization
methods, which may use various techniques. Third, the analysis of results is
possible via the post-processing of the computed time-dependent physical
quantities, and visualization of the flow properties. This is done via
implementing various packages and libraries that are standardized in the field
of computational astrophysics and supported by community developers. In the
present paper, we discuss the physical picture of the cosmic sources which are
modeled using numerical framework. We also describe several technical issues,
in the particular context of our own experience with the performance of the
GRMHD code which we develop. We also present a suite of performance tests, done
on the High-Performance Computer cluster (HPC) in the Center for Mathematical
Modeling of the Warsaw University.
