We showcase the port to Julia of a massively parallel Multi-GPU solver for spontaneous nonlinear multi-physics flow localization in 3-D. The original solver is itself the result of a translation from a Matlab prototype to CUDA C and MPI. Our contribution is an illustration of Julia solving "the two language problem": the Matlab prototype and the CUDA C + MPI production code are being replaced by a single Julia code that will serve both further prototyping and production. The solver's parallel and matrix-free design enables a short time to solution and is applicable to solve a wide variety of coupled and nonlinear systems of partial differential equations in 3-D. The employed stencil-based iterative method optimally suits both shared and distributed memory parallelization. As reference, the original Multi-GPU solver achieved a high performance and a nearly ideal parallel efficiency on up to 5120 NVIDIA Tesla P100 GPUs on the hybrid Cray XC-50 "Piz Daint" supercomputer at the Swiss National Supercomputing Centre, CSCS. We report the first performance and scaling results obtained with the Julia port. We present additionally our porting approach and discuss the related challenges.