 
     // copy back flux
     cudaMemcpy(F_h, F_d, nx*ny*vec_dim*sizeof(float), cudaMemcpyDeviceToHost);
-    //bcs_fv(F_h, nx, ny, nlayers, ng);
+
     if (n_processes == 1) {
         bcs_fv(F_h, nx, ny, ng, vec_dim);
     } else {


            Sy *= phi(r);

            //if (nan_check(q_comp[coarse_index])) printf("q_comp is nan :(\n");

            q_f[(2*y * nxf + 2*x) * 4 + n] = q_comp[coarse_index] -
                0.25 * (dx * Sx + dy * Sy);



    // copy back flux
    cudaMemcpy(F_h, F_d, nx*ny*vec_dim*sizeof(float), cudaMemcpyDeviceToHost);

    if (n_processes == 1) {
        bcs_fv(F_h, nx, ny, ng, vec_dim);
    } else {
