We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Compute initial streamfunction from initial velocity
Not knowing the numerics, and possibly not reading the docs carefully enough, it is unclear to me how to initialize the velocities in the model.  

For a channel run, re-entrant in x, with intial velocity 0.1 m/s everywhere, no forcing, I tried do in `set_initial_conditions`: `s.u = update(vs.u, at[...], 0.1 * vs.maskU[..., None])`.  

The velocity signal only lasts for one time step, and then it is gone.  It _does_ create a small pressure perturbations that drive internal waves, but the mean flow of 0.1 m/s is immediately gone.  Conversely, the initial conditions have psi=0 everywhere, and then immediately on the next time step there is a stream function, but if the units are really m^3/s it is far too small.  

Was I to initialize psi at the beginning instead of u, or in addition to u?  

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `veros/core/external/streamfunction_init.py`
Content:
```
1 from veros import logger, veros_kernel, veros_routine, KernelOutput
2 from veros.variables import allocate
3 from veros.distributed import global_max
4 from veros.core import utilities as mainutils
5 from veros.core.operators import numpy as npx, update, at
6 from veros.core.external import island, line_integrals
7 from veros.core.external.solvers import get_linear_solver
8 
9 
10 @veros_routine
11 def get_isleperim(state):
12     """
13     preprocess land map using MOMs algorithm for B-grid to determine number of islands
14     """
15     from veros.state import resize_dimension
16 
17     vs = state.variables
18 
19     island.isleperim(state)
20 
21     # now that we know the number of islands we can resize
22     # all arrays depending on that
23     nisle = int(global_max(npx.max(vs.land_map)))
24     resize_dimension(state, "isle", nisle)
25     vs.isle = npx.arange(nisle)
26 
27 
28 @veros_routine
29 def streamfunction_init(state):
30     """
31     prepare for island integrals
32     """
33     vs = state.variables
34     settings = state.settings
35 
36     logger.info("Initializing streamfunction method")
37 
38     get_isleperim(state)
39 
40     vs.update(boundary_masks(state))
41 
42     # populate linear solver cache
43     linear_solver = get_linear_solver(state)
44 
45     """
46     precalculate time independent boundary components of streamfunction
47     """
48     forc = allocate(state.dimensions, ("xt", "yt"))
49 
50     vs.psin = update(vs.psin, at[...], vs.maskZ[..., -1, npx.newaxis])
51 
52     for isle in range(state.dimensions["isle"]):
53         logger.info(f" Solving for boundary contribution by island {isle:d}")
54         isle_boundary = (
55             vs.line_dir_east_mask[..., isle]
56             | vs.line_dir_west_mask[..., isle]
57             | vs.line_dir_north_mask[..., isle]
58             | vs.line_dir_south_mask[..., isle]
59         )
60         isle_sol = linear_solver.solve(state, forc, vs.psin[:, :, isle], boundary_val=isle_boundary)
61         vs.psin = update(vs.psin, at[:, :, isle], isle_sol)
62 
63     vs.psin = mainutils.enforce_boundaries(vs.psin, settings.enable_cyclic_x)
64 
65     line_psin_out = island_integrals(state)
66     vs.update(line_psin_out)
67 
68 
69 @veros_kernel
70 def island_integrals(state):
71     """
72     precalculate time independent island integrals
73     """
74     vs = state.variables
75 
76     uloc = allocate(state.dimensions, ("xt", "yt", "isle"))
77     vloc = allocate(state.dimensions, ("xt", "yt", "isle"))
78 
79     uloc = update(
80         uloc,
81         at[1:, 1:, :],
82         -(vs.psin[1:, 1:, :] - vs.psin[1:, :-1, :])
83         * vs.maskU[1:, 1:, -1, npx.newaxis]
84         / vs.dyt[npx.newaxis, 1:, npx.newaxis]
85         * vs.hur[1:, 1:, npx.newaxis],
86     )
87 
88     vloc = update(
89         vloc,
90         at[1:, 1:, ...],
91         (vs.psin[1:, 1:, :] - vs.psin[:-1, 1:, :])
92         * vs.maskV[1:, 1:, -1, npx.newaxis]
93         / (vs.cosu[npx.newaxis, 1:, npx.newaxis] * vs.dxt[1:, npx.newaxis, npx.newaxis])
94         * vs.hvr[1:, 1:, npx.newaxis],
95     )
96 
97     vs.line_psin = line_integrals.line_integrals(state, uloc=uloc, vloc=vloc, kind="full")
98     return KernelOutput(line_psin=vs.line_psin)
99 
100 
101 @veros_kernel
102 def boundary_masks(state):
103     """
104     now that the number of islands is known we can allocate the rest of the variables
105     """
106     vs = state.variables
107     settings = state.settings
108 
109     boundary_map = vs.land_map[..., npx.newaxis] == npx.arange(1, state.dimensions["isle"] + 1)
110 
111     if settings.enable_cyclic_x:
112         vs.line_dir_east_mask = update(
113             vs.line_dir_east_mask, at[2:-2, 1:-1], boundary_map[3:-1, 1:-1] & ~boundary_map[3:-1, 2:]
114         )
115         vs.line_dir_west_mask = update(
116             vs.line_dir_west_mask, at[2:-2, 1:-1], boundary_map[2:-2, 2:] & ~boundary_map[2:-2, 1:-1]
117         )
118         vs.line_dir_south_mask = update(
119             vs.line_dir_south_mask, at[2:-2, 1:-1], boundary_map[2:-2, 1:-1] & ~boundary_map[3:-1, 1:-1]
120         )
121         vs.line_dir_north_mask = update(
122             vs.line_dir_north_mask, at[2:-2, 1:-1], boundary_map[3:-1, 2:] & ~boundary_map[2:-2, 2:]
123         )
124     else:
125         vs.line_dir_east_mask = update(
126             vs.line_dir_east_mask, at[1:-1, 1:-1], boundary_map[2:, 1:-1] & ~boundary_map[2:, 2:]
127         )
128         vs.line_dir_west_mask = update(
129             vs.line_dir_west_mask, at[1:-1, 1:-1], boundary_map[1:-1, 2:] & ~boundary_map[1:-1, 1:-1]
130         )
131         vs.line_dir_south_mask = update(
132             vs.line_dir_south_mask, at[1:-1, 1:-1], boundary_map[1:-1, 1:-1] & ~boundary_map[2:, 1:-1]
133         )
134         vs.line_dir_north_mask = update(
135             vs.line_dir_north_mask, at[1:-1, 1:-1], boundary_map[2:, 2:] & ~boundary_map[1:-1, 2:]
136         )
137 
138     vs.isle_boundary_mask = ~npx.any(
139         vs.line_dir_east_mask | vs.line_dir_west_mask | vs.line_dir_south_mask | vs.line_dir_north_mask, axis=2
140     )
141 
142     return KernelOutput(
143         isle_boundary_mask=vs.isle_boundary_mask,
144         line_dir_east_mask=vs.line_dir_east_mask,
145         line_dir_west_mask=vs.line_dir_west_mask,
146         line_dir_south_mask=vs.line_dir_south_mask,
147         line_dir_north_mask=vs.line_dir_north_mask,
148     )
149
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/veros/core/external/streamfunction_init.py b/veros/core/external/streamfunction_init.py
--- a/veros/core/external/streamfunction_init.py
+++ b/veros/core/external/streamfunction_init.py
@@ -3,7 +3,7 @@
 from veros.distributed import global_max
 from veros.core import utilities as mainutils
 from veros.core.operators import numpy as npx, update, at
-from veros.core.external import island, line_integrals
+from veros.core.external import island, line_integrals, solve_stream
 from veros.core.external.solvers import get_linear_solver
 
 
@@ -65,6 +65,29 @@
     line_psin_out = island_integrals(state)
     vs.update(line_psin_out)
 
+    """
+    take care of initial velocity
+    """
+
+    # transfer initial velocity to tendency
+    vs.du = update(vs.du, at[..., vs.tau], vs.u[..., vs.tau] / settings.dt_mom / (1.5 + settings.AB_eps))
+    vs.dv = update(vs.dv, at[..., vs.tau], vs.v[..., vs.tau] / settings.dt_mom / (1.5 + settings.AB_eps))
+    vs.u = update(vs.u, at[...], 0)
+    vs.v = update(vs.v, at[...], 0)
+
+    # run streamfunction solver to determine initial barotropic and baroclinic modes
+    solve_stream.solve_streamfunction(state)
+
+    vs.psi = update(vs.psi, at[...], vs.psi[..., vs.taup1, npx.newaxis])
+    vs.u = update(
+        vs.u, at[...], mainutils.enforce_boundaries(vs.u[..., vs.taup1, npx.newaxis], settings.enable_cyclic_x)
+    )
+    vs.v = update(
+        vs.v, at[...], mainutils.enforce_boundaries(vs.v[..., vs.taup1, npx.newaxis], settings.enable_cyclic_x)
+    )
+    vs.du = update(vs.du, at[..., vs.tau], 0)
+    vs.dv = update(vs.dv, at[..., vs.tau], 0)
+
 
 @veros_kernel
 def island_integrals(state):

{"golden_diff": "diff --git a/veros/core/external/streamfunction_init.py b/veros/core/external/streamfunction_init.py\n--- a/veros/core/external/streamfunction_init.py\n+++ b/veros/core/external/streamfunction_init.py\n@@ -3,7 +3,7 @@\n from veros.distributed import global_max\n from veros.core import utilities as mainutils\n from veros.core.operators import numpy as npx, update, at\n-from veros.core.external import island, line_integrals\n+from veros.core.external import island, line_integrals, solve_stream\n from veros.core.external.solvers import get_linear_solver\n \n \n@@ -65,6 +65,29 @@\n     line_psin_out = island_integrals(state)\n     vs.update(line_psin_out)\n \n+    \"\"\"\n+    take care of initial velocity\n+    \"\"\"\n+\n+    # transfer initial velocity to tendency\n+    vs.du = update(vs.du, at[..., vs.tau], vs.u[..., vs.tau] / settings.dt_mom / (1.5 + settings.AB_eps))\n+    vs.dv = update(vs.dv, at[..., vs.tau], vs.v[..., vs.tau] / settings.dt_mom / (1.5 + settings.AB_eps))\n+    vs.u = update(vs.u, at[...], 0)\n+    vs.v = update(vs.v, at[...], 0)\n+\n+    # run streamfunction solver to determine initial barotropic and baroclinic modes\n+    solve_stream.solve_streamfunction(state)\n+\n+    vs.psi = update(vs.psi, at[...], vs.psi[..., vs.taup1, npx.newaxis])\n+    vs.u = update(\n+        vs.u, at[...], mainutils.enforce_boundaries(vs.u[..., vs.taup1, npx.newaxis], settings.enable_cyclic_x)\n+    )\n+    vs.v = update(\n+        vs.v, at[...], mainutils.enforce_boundaries(vs.v[..., vs.taup1, npx.newaxis], settings.enable_cyclic_x)\n+    )\n+    vs.du = update(vs.du, at[..., vs.tau], 0)\n+    vs.dv = update(vs.dv, at[..., vs.tau], 0)\n+\n \n @veros_kernel\n def island_integrals(state):\n", "issue": "Compute initial streamfunction from initial velocity\nNot knowing the numerics, and possibly not reading the docs carefully enough, it is unclear to me how to initialize the velocities in the model.  \r\n\r\nFor a channel run, re-entrant in x, with intial velocity 0.1 m/s everywhere, no forcing, I tried do in `set_initial_conditions`: `s.u = update(vs.u, at[...], 0.1 * vs.maskU[..., None])`.  \r\n\r\nThe velocity signal only lasts for one time step, and then it is gone.  It _does_ create a small pressure perturbations that drive internal waves, but the mean flow of 0.1 m/s is immediately gone.  Conversely, the initial conditions have psi=0 everywhere, and then immediately on the next time step there is a stream function, but if the units are really m^3/s it is far too small.  \r\n\r\nWas I to initialize psi at the beginning instead of u, or in addition to u?  \n", "before_files": [{"content": "from veros import logger, veros_kernel, veros_routine, KernelOutput\nfrom veros.variables import allocate\nfrom veros.distributed import global_max\nfrom veros.core import utilities as mainutils\nfrom veros.core.operators import numpy as npx, update, at\nfrom veros.core.external import island, line_integrals\nfrom veros.core.external.solvers import get_linear_solver\n\n\n@veros_routine\ndef get_isleperim(state):\n    \"\"\"\n    preprocess land map using MOMs algorithm for B-grid to determine number of islands\n    \"\"\"\n    from veros.state import resize_dimension\n\n    vs = state.variables\n\n    island.isleperim(state)\n\n    # now that we know the number of islands we can resize\n    # all arrays depending on that\n    nisle = int(global_max(npx.max(vs.land_map)))\n    resize_dimension(state, \"isle\", nisle)\n    vs.isle = npx.arange(nisle)\n\n\n@veros_routine\ndef streamfunction_init(state):\n    \"\"\"\n    prepare for island integrals\n    \"\"\"\n    vs = state.variables\n    settings = state.settings\n\n    logger.info(\"Initializing streamfunction method\")\n\n    get_isleperim(state)\n\n    vs.update(boundary_masks(state))\n\n    # populate linear solver cache\n    linear_solver = get_linear_solver(state)\n\n    \"\"\"\n    precalculate time independent boundary components of streamfunction\n    \"\"\"\n    forc = allocate(state.dimensions, (\"xt\", \"yt\"))\n\n    vs.psin = update(vs.psin, at[...], vs.maskZ[..., -1, npx.newaxis])\n\n    for isle in range(state.dimensions[\"isle\"]):\n        logger.info(f\" Solving for boundary contribution by island {isle:d}\")\n        isle_boundary = (\n            vs.line_dir_east_mask[..., isle]\n            | vs.line_dir_west_mask[..., isle]\n            | vs.line_dir_north_mask[..., isle]\n            | vs.line_dir_south_mask[..., isle]\n        )\n        isle_sol = linear_solver.solve(state, forc, vs.psin[:, :, isle], boundary_val=isle_boundary)\n        vs.psin = update(vs.psin, at[:, :, isle], isle_sol)\n\n    vs.psin = mainutils.enforce_boundaries(vs.psin, settings.enable_cyclic_x)\n\n    line_psin_out = island_integrals(state)\n    vs.update(line_psin_out)\n\n\n@veros_kernel\ndef island_integrals(state):\n    \"\"\"\n    precalculate time independent island integrals\n    \"\"\"\n    vs = state.variables\n\n    uloc = allocate(state.dimensions, (\"xt\", \"yt\", \"isle\"))\n    vloc = allocate(state.dimensions, (\"xt\", \"yt\", \"isle\"))\n\n    uloc = update(\n        uloc,\n        at[1:, 1:, :],\n        -(vs.psin[1:, 1:, :] - vs.psin[1:, :-1, :])\n        * vs.maskU[1:, 1:, -1, npx.newaxis]\n        / vs.dyt[npx.newaxis, 1:, npx.newaxis]\n        * vs.hur[1:, 1:, npx.newaxis],\n    )\n\n    vloc = update(\n        vloc,\n        at[1:, 1:, ...],\n        (vs.psin[1:, 1:, :] - vs.psin[:-1, 1:, :])\n        * vs.maskV[1:, 1:, -1, npx.newaxis]\n        / (vs.cosu[npx.newaxis, 1:, npx.newaxis] * vs.dxt[1:, npx.newaxis, npx.newaxis])\n        * vs.hvr[1:, 1:, npx.newaxis],\n    )\n\n    vs.line_psin = line_integrals.line_integrals(state, uloc=uloc, vloc=vloc, kind=\"full\")\n    return KernelOutput(line_psin=vs.line_psin)\n\n\n@veros_kernel\ndef boundary_masks(state):\n    \"\"\"\n    now that the number of islands is known we can allocate the rest of the variables\n    \"\"\"\n    vs = state.variables\n    settings = state.settings\n\n    boundary_map = vs.land_map[..., npx.newaxis] == npx.arange(1, state.dimensions[\"isle\"] + 1)\n\n    if settings.enable_cyclic_x:\n        vs.line_dir_east_mask = update(\n            vs.line_dir_east_mask, at[2:-2, 1:-1], boundary_map[3:-1, 1:-1] & ~boundary_map[3:-1, 2:]\n        )\n        vs.line_dir_west_mask = update(\n            vs.line_dir_west_mask, at[2:-2, 1:-1], boundary_map[2:-2, 2:] & ~boundary_map[2:-2, 1:-1]\n        )\n        vs.line_dir_south_mask = update(\n            vs.line_dir_south_mask, at[2:-2, 1:-1], boundary_map[2:-2, 1:-1] & ~boundary_map[3:-1, 1:-1]\n        )\n        vs.line_dir_north_mask = update(\n            vs.line_dir_north_mask, at[2:-2, 1:-1], boundary_map[3:-1, 2:] & ~boundary_map[2:-2, 2:]\n        )\n    else:\n        vs.line_dir_east_mask = update(\n            vs.line_dir_east_mask, at[1:-1, 1:-1], boundary_map[2:, 1:-1] & ~boundary_map[2:, 2:]\n        )\n        vs.line_dir_west_mask = update(\n            vs.line_dir_west_mask, at[1:-1, 1:-1], boundary_map[1:-1, 2:] & ~boundary_map[1:-1, 1:-1]\n        )\n        vs.line_dir_south_mask = update(\n            vs.line_dir_south_mask, at[1:-1, 1:-1], boundary_map[1:-1, 1:-1] & ~boundary_map[2:, 1:-1]\n        )\n        vs.line_dir_north_mask = update(\n            vs.line_dir_north_mask, at[1:-1, 1:-1], boundary_map[2:, 2:] & ~boundary_map[1:-1, 2:]\n        )\n\n    vs.isle_boundary_mask = ~npx.any(\n        vs.line_dir_east_mask | vs.line_dir_west_mask | vs.line_dir_south_mask | vs.line_dir_north_mask, axis=2\n    )\n\n    return KernelOutput(\n        isle_boundary_mask=vs.isle_boundary_mask,\n        line_dir_east_mask=vs.line_dir_east_mask,\n        line_dir_west_mask=vs.line_dir_west_mask,\n        line_dir_south_mask=vs.line_dir_south_mask,\n        line_dir_north_mask=vs.line_dir_north_mask,\n    )\n", "path": "veros/core/external/streamfunction_init.py"}], "after_files": [{"content": "from veros import logger, veros_kernel, veros_routine, KernelOutput\nfrom veros.variables import allocate\nfrom veros.distributed import global_max\nfrom veros.core import utilities as mainutils\nfrom veros.core.operators import numpy as npx, update, at\nfrom veros.core.external import island, line_integrals, solve_stream\nfrom veros.core.external.solvers import get_linear_solver\n\n\n@veros_routine\ndef get_isleperim(state):\n    \"\"\"\n    preprocess land map using MOMs algorithm for B-grid to determine number of islands\n    \"\"\"\n    from veros.state import resize_dimension\n\n    vs = state.variables\n\n    island.isleperim(state)\n\n    # now that we know the number of islands we can resize\n    # all arrays depending on that\n    nisle = int(global_max(npx.max(vs.land_map)))\n    resize_dimension(state, \"isle\", nisle)\n    vs.isle = npx.arange(nisle)\n\n\n@veros_routine\ndef streamfunction_init(state):\n    \"\"\"\n    prepare for island integrals\n    \"\"\"\n    vs = state.variables\n    settings = state.settings\n\n    logger.info(\"Initializing streamfunction method\")\n\n    get_isleperim(state)\n\n    vs.update(boundary_masks(state))\n\n    # populate linear solver cache\n    linear_solver = get_linear_solver(state)\n\n    \"\"\"\n    precalculate time independent boundary components of streamfunction\n    \"\"\"\n    forc = allocate(state.dimensions, (\"xt\", \"yt\"))\n\n    vs.psin = update(vs.psin, at[...], vs.maskZ[..., -1, npx.newaxis])\n\n    for isle in range(state.dimensions[\"isle\"]):\n        logger.info(f\" Solving for boundary contribution by island {isle:d}\")\n        isle_boundary = (\n            vs.line_dir_east_mask[..., isle]\n            | vs.line_dir_west_mask[..., isle]\n            | vs.line_dir_north_mask[..., isle]\n            | vs.line_dir_south_mask[..., isle]\n        )\n        isle_sol = linear_solver.solve(state, forc, vs.psin[:, :, isle], boundary_val=isle_boundary)\n        vs.psin = update(vs.psin, at[:, :, isle], isle_sol)\n\n    vs.psin = mainutils.enforce_boundaries(vs.psin, settings.enable_cyclic_x)\n\n    line_psin_out = island_integrals(state)\n    vs.update(line_psin_out)\n\n    \"\"\"\n    take care of initial velocity\n    \"\"\"\n\n    # transfer initial velocity to tendency\n    vs.du = update(vs.du, at[..., vs.tau], vs.u[..., vs.tau] / settings.dt_mom / (1.5 + settings.AB_eps))\n    vs.dv = update(vs.dv, at[..., vs.tau], vs.v[..., vs.tau] / settings.dt_mom / (1.5 + settings.AB_eps))\n    vs.u = update(vs.u, at[...], 0)\n    vs.v = update(vs.v, at[...], 0)\n\n    # run streamfunction solver to determine initial barotropic and baroclinic modes\n    solve_stream.solve_streamfunction(state)\n\n    vs.psi = update(vs.psi, at[...], vs.psi[..., vs.taup1, npx.newaxis])\n    vs.u = update(\n        vs.u, at[...], mainutils.enforce_boundaries(vs.u[..., vs.taup1, npx.newaxis], settings.enable_cyclic_x)\n    )\n    vs.v = update(\n        vs.v, at[...], mainutils.enforce_boundaries(vs.v[..., vs.taup1, npx.newaxis], settings.enable_cyclic_x)\n    )\n    vs.du = update(vs.du, at[..., vs.tau], 0)\n    vs.dv = update(vs.dv, at[..., vs.tau], 0)\n\n\n@veros_kernel\ndef island_integrals(state):\n    \"\"\"\n    precalculate time independent island integrals\n    \"\"\"\n    vs = state.variables\n\n    uloc = allocate(state.dimensions, (\"xt\", \"yt\", \"isle\"))\n    vloc = allocate(state.dimensions, (\"xt\", \"yt\", \"isle\"))\n\n    uloc = update(\n        uloc,\n        at[1:, 1:, :],\n        -(vs.psin[1:, 1:, :] - vs.psin[1:, :-1, :])\n        * vs.maskU[1:, 1:, -1, npx.newaxis]\n        / vs.dyt[npx.newaxis, 1:, npx.newaxis]\n        * vs.hur[1:, 1:, npx.newaxis],\n    )\n\n    vloc = update(\n        vloc,\n        at[1:, 1:, ...],\n        (vs.psin[1:, 1:, :] - vs.psin[:-1, 1:, :])\n        * vs.maskV[1:, 1:, -1, npx.newaxis]\n        / (vs.cosu[npx.newaxis, 1:, npx.newaxis] * vs.dxt[1:, npx.newaxis, npx.newaxis])\n        * vs.hvr[1:, 1:, npx.newaxis],\n    )\n\n    vs.line_psin = line_integrals.line_integrals(state, uloc=uloc, vloc=vloc, kind=\"full\")\n    return KernelOutput(line_psin=vs.line_psin)\n\n\n@veros_kernel\ndef boundary_masks(state):\n    \"\"\"\n    now that the number of islands is known we can allocate the rest of the variables\n    \"\"\"\n    vs = state.variables\n    settings = state.settings\n\n    boundary_map = vs.land_map[..., npx.newaxis] == npx.arange(1, state.dimensions[\"isle\"] + 1)\n\n    if settings.enable_cyclic_x:\n        vs.line_dir_east_mask = update(\n            vs.line_dir_east_mask, at[2:-2, 1:-1], boundary_map[3:-1, 1:-1] & ~boundary_map[3:-1, 2:]\n        )\n        vs.line_dir_west_mask = update(\n            vs.line_dir_west_mask, at[2:-2, 1:-1], boundary_map[2:-2, 2:] & ~boundary_map[2:-2, 1:-1]\n        )\n        vs.line_dir_south_mask = update(\n            vs.line_dir_south_mask, at[2:-2, 1:-1], boundary_map[2:-2, 1:-1] & ~boundary_map[3:-1, 1:-1]\n        )\n        vs.line_dir_north_mask = update(\n            vs.line_dir_north_mask, at[2:-2, 1:-1], boundary_map[3:-1, 2:] & ~boundary_map[2:-2, 2:]\n        )\n    else:\n        vs.line_dir_east_mask = update(\n            vs.line_dir_east_mask, at[1:-1, 1:-1], boundary_map[2:, 1:-1] & ~boundary_map[2:, 2:]\n        )\n        vs.line_dir_west_mask = update(\n            vs.line_dir_west_mask, at[1:-1, 1:-1], boundary_map[1:-1, 2:] & ~boundary_map[1:-1, 1:-1]\n        )\n        vs.line_dir_south_mask = update(\n            vs.line_dir_south_mask, at[1:-1, 1:-1], boundary_map[1:-1, 1:-1] & ~boundary_map[2:, 1:-1]\n        )\n        vs.line_dir_north_mask = update(\n            vs.line_dir_north_mask, at[1:-1, 1:-1], boundary_map[2:, 2:] & ~boundary_map[1:-1, 2:]\n        )\n\n    vs.isle_boundary_mask = ~npx.any(\n        vs.line_dir_east_mask | vs.line_dir_west_mask | vs.line_dir_south_mask | vs.line_dir_north_mask, axis=2\n    )\n\n    return KernelOutput(\n        isle_boundary_mask=vs.isle_boundary_mask,\n        line_dir_east_mask=vs.line_dir_east_mask,\n        line_dir_west_mask=vs.line_dir_west_mask,\n        line_dir_south_mask=vs.line_dir_south_mask,\n        line_dir_north_mask=vs.line_dir_north_mask,\n    )\n", "path": "veros/core/external/streamfunction_init.py"}]}