Test river exchange XMI

[HendrikKok]

#%%
import os
import datetime
import numpy as np
from xmipy import XmiWrapper

os.chdir(r"c:\werkmap\test_mf6_xmi")

# Type hints
FloatArray = np.ndarray
IntArray = np.ndarray
BoolArray = np.ndarray

#%%



class Simulation:
    """
    Run all stress periods in a simulation 
    """
    def __init__(self, wdir: str, name: str):
        self.modelname = name
        self.mf6 = XmiWrapper(lib_path="libmf6.dll", working_directory=wdir)
        self.mf6.initialize()
        self.max_iter = self.mf6.get_value_ptr("SLN_1/MXITER")[0]
        shape = np.zeros(1, dtype=np.int32)
        self.ncell = self.mf6.get_grid_shape(1, shape)[0]
        #(f"Initialized model with {self.ncell} cells")

    def do_iter(self, sol_id: int) -> bool:
        """Execute a single iteration"""
        has_converged = self.mf6.solve(sol_id)
        return has_converged

    def update(self):
        # We cannot set the timestep (yet) in Modflow
        # -> set to the (dummy) value 0.0 for now
        self.mf6.prepare_time_step(0.0)
        #print(self.mf6.get_value_ptr("MODEL/RIV-1/BOUND"),'RIV-bound')
        #print(self.mf6.get_value_ptr("MODEL/RIV-1/NODELIST"),'RIV-MFid')
        #print(self.mf6.get_value_ptr("SLN_1/X"),'head')
        self.mf6.prepare_solve(1)
        # Convergence loop
        for kiter in range(1, self.max_iter + 1):
            #(f"MF6 outer iteration: {kiter}")
            has_converged = self.do_iter(1)
            if has_converged:
                #(f"MF6 converged in {kiter} iterations")
                break
        self.mf6.finalize_solve(1)
        # Finish timestep
        self.mf6.finalize_time_step()
        current_time = self.mf6.get_current_time()
        riv=self.mf6.get_value_ptr("MODEL/RIV-1/BOUND")
        riv_id=self.mf6.get_value_ptr("MODEL/RIV-1/NODELIST")
        head=self.mf6.get_value_ptr("SLN_1/X")
        iriv=-1
        q=0
        for inode in riv_id:
            iriv=+1
            q=q+(riv[iriv,1]*(riv[iriv,0]-head[inode]))
            
            
        print(current_time,np.sum(q))        
            
            
        if(current_time==10.0):
            riv[:,0]=np.array([-0.5,-0.5,-0.5])
            self.mf6.set_value("MODEL/RIV-1/BOUND",riv)   
            
        return current_time

    def get_times(self):
        """Return times"""
        return (
            self.mf6.get_start_time(),
            self.mf6.get_current_time(),
            self.mf6.get_end_time(),
        )

    def run(self, periods):
        iperiod = 0
        _, current_time, end_time = self.get_times()
        while (current_time < end_time) and iperiod < periods:
            #print(f"MF6 starting period {iperiod}")
            current_time = self.update()
            iperiod += 1
        print(f"Simulation terminated normally for {periods} periods")

    def finalize(self):
        self.mf6.finalize()



def run_3d_model(periods):
    wdir = "model"
    name = "GWF_1"
    sim = Simulation(wdir, name)
    start = datetime.datetime.now()
    sim.run(periods)
    end = datetime.datetime.now()
    print(end - start)
    sim.finalize()


run_3d_model(20)