Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

does Mf6Splitter work for multilayer models? #2415

Open
ougx opened this issue Jan 12, 2025 · 5 comments · May be fixed by #2418
Open

does Mf6Splitter work for multilayer models? #2415

ougx opened this issue Jan 12, 2025 · 5 comments · May be fixed by #2418

Comments

@ougx
Copy link
Contributor

ougx commented Jan 12, 2025

it does not seem to work for me with the version 3.9.1.

Below is a test:

import flopy
import os
import numpy as np
#%% create model
origname = "Orig"
os.makedirs(origname, exist_ok=True)
sim = flopy.mf6.MFSimulation(sim_ws=origname)


tdis = flopy.mf6.ModflowTdis(
    sim, nper=1, perioddata=[(1.0, 1, 1.0)]
)

ims = flopy.mf6.ModflowIms(sim)

gwf = flopy.mf6.ModflowGwf(
    sim,
    modelname=origname,
)

dis = flopy.mf6.ModflowGwedis(
    gwf,
    nlay=2,
    nrow=1,
    ncol=2,
    delr=1,
    delc=1,
    top=2,
    botm=[1,0]
)

ic = flopy.mf6.ModflowGwfic(gwf, strt=2.0)

npf = flopy.mf6.ModflowGwfnpf(
    gwf,
    k=1.0,
    k33=0.1,
)

wel = flopy.mf6.ModflowGwfwel(
    gwf,
    print_input=True,
    print_flows=True,
    stress_period_data={0:[[(0,0,0),0.001], [(0,0,1),0.001]]},
    save_flows=False,
)

sim.write_simulation()

#%% parallel
array = np.ones([2, 1, 2], int)
array[:,0,1] = 2

mfsplit = flopy.mf6.utils.Mf6Splitter(sim)
new_sim = mfsplit.split_model(array)
new_sim.set_sim_path(".")
new_sim.write_simulation()

#%% parallel one layer mask
array = np.ones([1, 2], int)
array[0,1] = 2

mfsplit = flopy.mf6.utils.Mf6Splitter(sim)
new_sim = mfsplit.split_model(array)
new_sim.set_sim_path(".")
new_sim.write_simulation()
@jlarsen-usgs
Copy link
Contributor

@ougx
The model splitter supports multi-layer models. It does not support splitting models differently by layer (e.g., splitting mask is different for layer 1 and layer 2). As such, the input array for splitting a model should have exactly (NROW, NCOL) or (NCPL) entries for a structured or vertex grid model, whether it's one layered or multi-layered. For an unstructured grid model the dimension of the splitting array should be NNODES.

For your model to work you would need to change:

array = np.ones([2, 1, 2], int)
array[:,0,1] = 2

to

array = np.ones([1, 1, 2], int)
array[:,0,1] = 2

or

array = np.ones([1, 2], int)
array[0, 1] = 2

For better load balancing with a large model, I'd recommend using the option:

mfsplit = flopy.mf6.utils.Mf6Splitter(sim)
mask = mfsplit.optimize_splitting_mask(2)
new_sim = mfsplit.split_model(mask)

@ougx
Copy link
Contributor Author

ougx commented Jan 13, 2025

Hi @jlarsen-usgs: thanks for your answer. But it doesn't seem to work for me:

array = np.ones([1, 1, 2], int)
array[:,0,1] = 2

mfsplit = flopy.mf6.utils.Mf6Splitter(sim)
mask = mfsplit.optimize_splitting_mask(2)
new_sim = mfsplit.split_model(array)
new_sim.set_sim_path(".")
new_sim.write_simulation()


Traceback (most recent call last):

  Cell In[387], line 6
    new_sim = mfsplit.split_model(array)

  File C:\Miniconda3\lib\site-packages\flopy\mf6\utils\model_splitter.py:3435 in split_model
    paks = self._remap_package(package)

  File C:\Miniconda3\lib\site-packages\flopy\mf6\utils\model_splitter.py:2971 in _remap_package
    mapped_data = self._remap_array(item, value, mapped_data)

  File C:\Miniconda3\lib\site-packages\flopy\mf6\utils\model_splitter.py:1204 in _remap_array
    new_array[new_nodes] = original_arr[old_nodes]

IndexError: index 1 is out of bounds for axis 0 with size 1

or

array = np.ones([1, 2], int)
array[0,1] = 2

mfsplit = flopy.mf6.utils.Mf6Splitter(sim)
mask = mfsplit.optimize_splitting_mask(2)
new_sim = mfsplit.split_model(array)
new_sim.set_sim_path(".")
new_sim.write_simulation()


Traceback (most recent call last):

  Cell In[388], line 6
    new_sim = mfsplit.split_model(array)

  File C:\Miniconda3\lib\site-packages\flopy\mf6\utils\model_splitter.py:3435 in split_model
    paks = self._remap_package(package)

  File C:\Miniconda3\lib\site-packages\flopy\mf6\utils\model_splitter.py:2971 in _remap_package
    mapped_data = self._remap_array(item, value, mapped_data)

  File C:\Miniconda3\lib\site-packages\flopy\mf6\utils\model_splitter.py:1204 in _remap_array
    new_array[new_nodes] = original_arr[old_nodes]

IndexError: index 1 is out of bounds for axis 0 with size 1

@jlarsen-usgs
Copy link
Contributor

@ougx

I see the issue now. I haven't wired in all of the support for GWE models yet, so the code itself it not recognizing that the GWE DIS package is a DIS package. Currently GWF and GWT are fully supported. I'll try to make the updates to support GWE models later today or tomorrow morning.

@jlarsen-usgs
Copy link
Contributor

I just noticed that you are building a GWF model and not a GWE model. As such, you should change your dis object to:

dis = flopy.mf6.ModflowGwfdis(
    gwf,
    nlay=2,
    nrow=1,
    ncol=2,
    delr=1,
    delc=1,
    top=2,
    botm=[1,0]
)

With that and using the code changes I had previously suggested, the splitter will work. Once the open PR is merged GWE packages will also be supported.

@ougx
Copy link
Contributor Author

ougx commented Jan 20, 2025

Ahh, my bad, Good catch. @jlarsen-usgs!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants