perfusion_micro.py

r"""
Homogenization of the Darcy flow in a thin porous layer.
The reference cell is composed of the matrix representing the dual porosity
and of two disconnected channels representing the primary porosity,
see paper [1].

[1] http://seth.asc.tuwien.ac.at/proc12/full_paper/Contribution183.pdf
"""

from __future__ import absolute_import
from sfepy.discrete.fem.periodic import match_x_plane, match_y_plane
import sfepy.homogenization.coefs_base as cb
import numpy as nm
from sfepy import data_dir
import six
from six.moves import range

def get_mats(pk, ph, pe, dim):

    m1 = nm.eye(dim, dtype=nm.float64) * pk
    m1[-1,-1] = pk / ph
    m2 = nm.eye(dim, dtype=nm.float64) * pk
    m2[-1,-1] = pk / ph ** 2

    return m1, m2

def recovery_perf(pb, corrs, macro):

    from sfepy.homogenization.recovery import compute_p_from_macro
    from sfepy.base.base import Struct

    slev = ''

    micro_coors = pb.domain.mesh.coors
    micro_nnod = pb.domain.mesh.n_nod

    centre_Y = nm.sum(pb.domain.mesh.coors, axis=0) / micro_nnod
    nodes_Y = {}

    channels = {}
    for k in six.iterkeys(macro):
        if 'press' in k:
            channels[k[-1]] = 1

    channels = list(channels.keys())

    varnames = ['pM']
    for ch in channels:
        nodes_Y[ch] = pb.domain.regions['Y' + ch].vertices
        varnames.append('p' + ch)

    pvars = pb.create_variables(varnames)

    press = {}

    # matrix
    press['M'] = \
       corrs['corrs_%s_gamma_p' % pb_def['name']]['pM'] * macro['g_p'] + \
       corrs['corrs_%s_gamma_m' % pb_def['name']]['pM'] * macro['g_m']

    out = {}
    # channels
    for ch in channels:
        press_mac = macro['press' + ch][0,0]
        press_mac_grad = macro['pressg' + ch]
        nnod = corrs['corrs_%s_pi%s' % (pb_def['name'], ch)]\
          ['p%s_0' % ch].shape[0]

        press_mic = nm.zeros( (nnod,1) )
        for key, val in \
          six.iteritems(corrs['corrs_%s_pi%s' % (pb_def['name'], ch)]):
            kk = int( key[-1] )
            press_mic += val * press_mac_grad[kk,0]

        for key in six.iterkeys(corrs):
            if ('_gamma_' + ch in key):
                kk = int(key[-1]) - 1
                press_mic += corrs[key]['p' + ch] * macro['g' + ch][kk]

        press_mic += \
          compute_p_from_macro(press_mac_grad[nm.newaxis,nm.newaxis,:,:],
                               micro_coors[nodes_Y[ch]], 0,
                               centre=centre_Y, extdim=-1).reshape((nnod,1))

        press[ch] = press_mac + eps0 * press_mic

        out[slev + 'p' + ch] = Struct(name='output_data',
                                      mode='vertex',
                                      data=press[ch],
                                      var_name='p' + ch,
                                      dofs=None)

        pvars['p' + ch].set_data(press_mic)
        dvel = pb.evaluate('ev_diffusion_velocity.iV.Y%s(mat1%s.k, p%s)'\
                           % (ch, ch, ch),
                           var_dict = {'p' + ch: pvars['p' + ch]},
                           mode='el_avg')

        out[slev + 'w' + ch] = Struct(name='output_data',
                                      mode='cell',
                                      data=dvel,
                                      var_name='w' + ch,
                                      dofs=None)

        press['M'] += corrs['corrs_%s_eta%s' % (pb_def['name'], ch)]['pM']\
          * press_mac

    pvars['pM'].set_data(press['M'])
    dvel = pb.evaluate('%e * ev_diffusion_velocity.iV.YM(mat1M.k, pM)' % eps0,
                       var_dict = {'pM': pvars['pM']}, mode='el_avg')

    out[slev + 'pM'] = Struct(name='output_data',
                              mode='vertex',
                              dat =press['M'],
                              var_name='pM',
                              dofs=None)

    out[slev + 'wM'] = Struct(name='output_data',
                              mode='cell',
                              data=dvel,
                              var_name='wM',
                              dofs=None )

    return out

geoms = {
    '2_4': ['2_4_Q1', '2', 5],
    '3_8': ['3_8_Q1', '4', 5],
    '3_4': ['3_4_P1', '3', 3],
    }

pb_def = {
    'name': '3d_2ch',
    'mesh_filename': data_dir + '/meshes/3d/perfusion_micro3d.mesh',
    'dim': 3,
    'geom': geoms['3_4'],
    'eps0': 1.0e-2,
    'param_h': 1.0,
    'param_kappa_m': 0.1,
    'matrix_mat_el_grp': 3,
    'channels': {
        'A': {
            'mat_el_grp': 1,
            'fix_nd_grp': (4, 1),
            'io_nd_grp': [ 1, 2, 3 ],
            'param_kappa_ch': 1.0,
            },
        'B': {
            'mat_el_grp': 2,
            'fix_nd_grp': (14, 11),
            'io_nd_grp': [ 11, 12, 13 ],
            'param_kappa_ch': 2.0,
            },
    },
}

filename_mesh = pb_def['mesh_filename']
eps0 = pb_def['eps0']
param_h = pb_def['param_h']

# integrals
integrals = {
    'iV' : 2,
    'iS' : 2,
}

functions = {
    'match_x_plane': (match_x_plane,),
    'match_y_plane': (match_y_plane,),
    }

aux = []
for ch, val in six.iteritems(pb_def['channels']):
    aux.append( 'r.bYM' + ch )

# basic regions
regions = {
    'Y': 'all',
    'YM': 'cells of group %d' % pb_def['matrix_mat_el_grp'],
    # periodic boundaries
    'Pl': ('vertices in (x < 0.001)', 'facet'),
    'Pr': ('vertices in (x > 0.999)', 'facet'),
    'PlYM': ('r.Pl *v r.YM', 'facet'),
    'PrYM': ('r.Pr *v r.YM', 'facet'),
    'bYMp': ('r.bYp *v r.YM', 'facet', 'YM'),
    'bYMm': ('r.bYm *v r.YM', 'facet', 'YM'),
    'bYMpm': ('r.bYMp +v r.bYMm', 'facet', 'YM'),
}

# matrix/channel boundaries
regions.update({
    'bYMchs': (' +v '.join(aux), 'facet', 'YM'),
    'YMmchs': 'r.YM -v r.bYMchs',
    })

# boundary conditions Gamma+/-
ebcs = {
    'gamma_pm_bYMchs': ('bYMchs', {'pM.0': 0.0}),
    'gamma_pm_YMmchs': ('YMmchs', {'pM.0': 1.0}),
    }

# periodic boundary conditions - matrix, X-direction
epbcs = {'periodic_xYM': (['PlYM', 'PrYM'], {'pM.0': 'pM.0'}, 'match_x_plane')}
lcbcs = {}

all_periodicYM = ['periodic_%sYM' % ii for ii in ['x', 'y'][:pb_def['dim']-1] ]
all_periodicY = {}

if pb_def['dim'] == 2:
    regions.update( {
        'bYm': ('vertices in (y < 0.001)', 'facet'),
        'bYp':  ('vertices in (y > 0.999)', 'facet'),
        } )
if pb_def['dim'] == 3:
    regions.update( {
        'Pn': ('vertices in (y < 0.001)', 'facet'),
        'Pf': ('vertices in (y > 0.999)', 'facet'),
        'PnYM': ('r.Pn *v r.YM', 'facet'),
        'PfYM': ('r.Pf *v r.YM', 'facet'),
        'bYm': ('vertices in (z < 0.001)', 'facet'),
        'bYp':  ('vertices in (z > 0.999)', 'facet'),
        } )
    # periodic boundary conditions - matrix, Y-direction
    epbcs.update( {
        'periodic_yYM': (['PnYM', 'PfYM'], {'pM.0': 'pM.0'}, 'match_y_plane'),
        } )

reg_io = {}
ebcs_eta = {}
ebcs_gamma = {}

# generate regions, ebcs, epbcs
for ch, val in six.iteritems(pb_def['channels']):

    all_periodicY[ch] = ['periodic_%sY%s' % (ii, ch)\
                         for ii in ['x', 'y'][:pb_def['dim']-1] ]

    # channels: YA, fixedYA, bYMA (matrix/channel boundaries)
    regions.update( {
        'Y' + ch: 'cells of group %d' % val['mat_el_grp'],
        'bYM' + ch: ('r.YM *v r.Y' + ch, 'facet', 'YM'),
        'PlY' + ch: ('r.Pl *v r.Y' + ch, 'facet'),
        'PrY' + ch: ('r.Pr *v r.Y' + ch, 'facet'),
        } )

    if 'fix_nd_grp' in val:
        regions.update({
            'fixedY' + ch: ('vertices of group %d' % val['fix_nd_grp'][0],
                            'vertex'),
            })

    ebcs_eta[ch] = []
    for ch2, val2 in six.iteritems(pb_def['channels']):
        aux = 'eta%s_bYM%s' % (ch, ch2)
        if ch2 == ch:
            ebcs.update( {aux: ('bYM' + ch2, {'pM.0': 1.0})} )
        else:
            ebcs.update( {aux: ('bYM' + ch2, {'pM.0': 0.0})} )

        ebcs_eta[ch].append(aux)

    # boundary conditions
    # periodic boundary conditions - channels, X-direction
    epbcs.update({
            'periodic_xY' + ch: (['PlY' + ch, 'PrY' + ch],
                                 {'p%s.0' % ch: 'p%s.0' % ch},
                                 'match_x_plane'),
            })

    if pb_def['dim'] == 3:
        regions.update({
                'PnY' + ch: ('r.Pn *v r.Y' + ch, 'facet'),
                'PfY' + ch: ('r.Pf *v r.Y' + ch, 'facet'),
                })
        # periodic boundary conditions - channels, Y-direction
        epbcs.update({
                'periodic_yY' + ch: (['PnY' + ch, 'PfY' + ch],
                                     {'p%s.0' % ch: 'p%s.0' % ch},
                                     'match_y_plane'),
            })

    reg_io[ch] = []
    aux_bY = []
    # channel: inputs/outputs
    for i_io in range( len( val['io_nd_grp'] ) ):
        io = '%s_%d' % (ch, i_io+1)

        # regions
        aux = val['io_nd_grp'][i_io]
        if 'fix_nd_grp' in val and val['fix_nd_grp'][1] == aux:
            regions.update( {
                'bY%s' % io : ('vertices of group %d +v r.fixedY%s' % (aux, ch),
                               'facet', 'Y%s' % ch),
                } )
        else:
            regions.update( {
                'bY%s' % io : ('vertices of group %d' % aux,
                               'facet', 'Y%s' % ch),
            } )

        aux_bY.append('r.bY%s' % io)
        reg_io[ch].append('bY%s' % io)

    regions.update({
        'bY' + ch : (' +v '.join(aux_bY), 'facet', 'Y' + ch),
        })

    # channel: inputs/outputs
    for i_io in range( len( val['io_nd_grp'] ) ):
        io = '%s_%d' % (ch, i_io + 1)
        ion = '%s_n%d' % (ch, i_io + 1)
        regions.update({
            'bY%s' % ion: ('r.bY%s -v r.bY%s' % (ch, io), 'facet', 'Y%s' % ch),
            })

        # boundary conditions
        aux = 'fix_p%s_bY%s' % (ch, ion)
        ebcs.update({
            aux: ('bY%s' % ion, {'p%s.0' % ch: 0.0}),
            })

    lcbcs.update({
        'imv' + ch: ('Y' + ch, {'ls%s.all' % ch: None}, None,
                     'integral_mean_value')
        })

###########################################################################
# materials, fields, variables, integrals

matk1, matk2 = get_mats(pb_def['param_kappa_m'], param_h, eps0, pb_def['dim'])

materials = {
    'mat1M': ({'k': matk1},),
    'mat2M': ({'k': matk2},),
}

fields = {
    'corrector_M': ('real', 'scalar', 'YM', 1),
    'vel_M': ('real', 'vector', 'YM', 1),
    'vol_all': ('real', 'scalar', 'Y', 1),
}

variables = {
    'pM': ('unknown field', 'corrector_M'),
    'qM': ('test field', 'corrector_M', 'pM'),
    'Pi_M': ('parameter field', 'corrector_M', '(set-to-None)'),
    'corr_M': ('parameter field', 'corrector_M', '(set-to-None)'),
    'corr1_M': ('parameter field', 'corrector_M', '(set-to-None)'),
    'corr2_M': ('parameter field', 'corrector_M', '(set-to-None)'),
    'wM': ('parameter field', 'vel_M', '(set-to-None)'),
    'vol_all': ('parameter field', 'vol_all', '(set-to-None)'),
}

# generate regions for channel inputs/outputs
for ch, val in six.iteritems(pb_def['channels']):

    matk1, matk2 = get_mats(val['param_kappa_ch'],  param_h,
                            eps0, pb_def['dim'])
    materials.update({
        'mat1' + ch: ({'k': matk1},),
        'mat2' + ch: ({'k': matk2},),
        })

    fields.update({
        'corrector_' + ch: ('real', 'scalar', 'Y' + ch, 1),
        'vel_' + ch: ('real', 'vector', 'Y' + ch, 1),
        })

    variables.update({
        'p' + ch: ('unknown field', 'corrector_' + ch),
        'q' + ch: ('test field', 'corrector_' + ch, 'p' + ch),
        'Pi_' + ch: ('parameter field', 'corrector_' + ch, '(set-to-None)'),
        'corr1_' + ch: ('parameter field', 'corrector_' + ch, '(set-to-None)'),
        'corr2_' + ch: ('parameter field', 'corrector_' + ch, '(set-to-None)'),
        'w' + ch: ('unknown field', 'vel_' + ch),
        # lagrange mutltipliers - integral mean value
        'ls' + ch: ('unknown field', 'corrector_' + ch),
        'lv' + ch: ('test field', 'corrector_' + ch, 'ls' + ch),
        })

###########################################################################

options = {
    'coefs': 'coefs',
    'requirements': 'requirements',
    'ls': 'ls', # linear solver to use
    'volumes': {
        'total': {
            'variables': ['vol_all'],
            'expression': """d_volume.iV.Y(vol_all)""",
            },
        'one': {
            'value': 1.0,
            }
        },
    'output_dir': './output',
    'file_per_var': True,
    'coefs_filename': 'coefs_perf_' + pb_def['name'],
    'coefs_info': {'eps0': eps0},
    'recovery_hook': 'recovery_perf',
    }

for ipm in ['p', 'm']:
    options['volumes'].update({
        'bYM' + ipm: {
            'variables': ['pM'],
            'expression': "d_surface.iS.bYM%s(pM)" % ipm,
            },
        'bY' + ipm: {
            'variables': ['vol_all'],
            'expression': "d_surface.iS.bY%s(vol_all)" % ipm,
            }
        })

for ch in six.iterkeys(reg_io):
    for ireg in reg_io[ch]:
        options['volumes'].update({
            ireg: {
                'variables': ['p' + ch],
                'expression': "d_surface.iS.%s(p%s)" % (ireg, ch),
                }
            })

###########################################################################
# equations, correctors, coefficients

coefs = {
    'vol_bYMpm': {
        'regions': ['bYMp', 'bYMm'],
        'expression': 'd_surface.iS.%s(pM)',
        'class': cb.VolumeFractions,
        },
    'filenames': {},
    }

# requirements for perfusion homog. coefficients
requirements = {
    'corrs_one_YM': {
        'variable': ['pM'],
        'ebcs': ['gamma_pm_YMmchs', 'gamma_pm_bYMchs'],
        'epbcs': [],
        'save_name': 'corrs_one_YM',
        'class': cb.CorrSetBCS,
        'dump_variables': ['pM'],
        },
    }

for ipm in ['p', 'm']:
    requirements.update({
        'corrs_gamma_' + ipm: {
            'requires': [],
            'ebcs': ['gamma_pm_bYMchs'],
            'epbcs': all_periodicYM,
            'equations': {
                'eq_gamma_pm': """dw_diffusion.iV.YM(mat2M.k, qM, pM) =
                             %e * dw_surface_integrate.iS.bYM%s(qM)"""\
                             % (1.0/param_h, ipm),
                },
            'class': cb.CorrOne,
            'save_name': 'corrs_%s_gamma_%s' % (pb_def['name'], ipm),
            'dump_variables': ['pM'],
            },
        })

    for ipm2 in ['p', 'm']:
        coefs.update({
            'H' + ipm + ipm2: { # test+
                'requires': ['corrs_gamma_' + ipm],
                'set_variables': [('corr_M', 'corrs_gamma_' + ipm, 'pM')],
                'expression': 'ev_surface_integrate.iS.bYM%s(corr_M)' % ipm2,
                'set_volume': 'bYp',
                'class': cb.CoefOne,
                },
        })

def get_channel(keys, bn):
    for ii in keys:
        if bn in ii:
            return ii[(ii.rfind(bn) + len(bn)):]

    return None

def set_corrpis(variables, ir, ic, mode, **kwargs):
    ch = get_channel(list(kwargs.keys()), 'pis_')
    pis = kwargs['pis_' + ch]
    corrs_pi = kwargs['corrs_pi' + ch]

    if mode == 'row':
        val = pis.states[ir]['p' + ch] + corrs_pi.states[ir]['p' + ch]
        variables['corr1_' + ch].set_data(val)
    elif mode == 'col':
        val = pis.states[ic]['p' + ch] + corrs_pi.states[ic]['p' + ch]
        variables['corr2_' + ch].set_data(val)

def set_corr_S(variables, ir, **kwargs):
    ch = get_channel(list(kwargs.keys()), 'pis_')
    io = get_channel(list(kwargs.keys()), 'corrs_gamma_')

    pis = kwargs['pis_' + ch]
    corrs_gamma = kwargs['corrs_gamma_' + io]

    pi = pis.states[ir]['p' + ch]
    val = corrs_gamma.state['p' + ch]
    variables['corr1_' + ch].set_data(pi)
    variables['corr2_' + ch].set_data(val)

def set_corr_cc(variables, ir, **kwargs):
    ch = get_channel(list(kwargs.keys()), 'pis_')
    pis = kwargs['pis_' + ch]
    corrs_pi = kwargs['corrs_pi' + ch]

    pi = pis.states[ir]['p' + ch]
    pi = pi - nm.mean(pi)
    val = pi + corrs_pi.states[ir]['p' + ch]
    variables['corr1_' + ch].set_data(val)

for ch, val in six.iteritems(pb_def['channels']):

    coefs.update({
        'G' + ch: {  # test+
            'requires': ['corrs_one' + ch, 'corrs_eta' + ch],
            'set_variables': [('corr1_M', 'corrs_one' + ch, 'pM'),
                              ('corr2_M', 'corrs_eta' + ch, 'pM')],
            'expression': 'dw_diffusion.iV.YM(mat2M.k, corr1_M, corr2_M)',
            'class': cb.CoefOne,
            },
        'K' + ch: {  # test+
            'requires': ['pis_' + ch, 'corrs_pi' + ch],
            'set_variables': set_corrpis,
            'expression': 'dw_diffusion.iV.Y%s(mat2%s.k, corr1_%s, corr2_%s)'\
                           % ((ch,) * 4),
            'dim': pb_def['dim'] - 1,
            'class': cb.CoefDimDim,
            },
        })

    requirements.update({
        'pis_' + ch: {
            'variables': ['p' + ch],
            'class': cb.ShapeDim,
            },
        'corrs_one' + ch: {
            'variable': ['pM'],
            'ebcs': ebcs_eta[ch],
            'epbcs': [],
            'save_name': 'corrs_%s_one%s' % (pb_def['name'], ch),
            'dump_variables': ['pM'],
            'class': cb.CorrSetBCS,
             },
        'corrs_eta' + ch: {
            'ebcs': ebcs_eta[ch],
            'epbcs': all_periodicYM,
            'equations': {
                'eq_eta': 'dw_diffusion.iV.YM(mat2M.k, qM, pM) = 0',
                },
            'class': cb.CorrOne,
            'save_name': 'corrs_%s_eta%s' % (pb_def['name'], ch),
            'dump_variables': ['pM'],
            },
        'corrs_pi' + ch: {
            'requires': ['pis_' + ch],
            'set_variables': [('Pi_' + ch, 'pis_' + ch, 'p' + ch)],
            'ebcs': [],
            'epbcs': all_periodicY[ch],
            'lcbcs': ['imv' + ch],
            'equations': {
                'eq_pi': """dw_diffusion.iV.Y%s(mat2%s.k, q%s, p%s)
                            + dw_volume_dot.iV.Y%s(q%s, ls%s)
                            = - dw_diffusion.iV.Y%s(mat2%s.k, q%s, Pi_%s)"""\
                            % ((ch,) * 11),
                'eq_imv': 'dw_volume_dot.iV.Y%s(lv%s, p%s) = 0' % ((ch,) * 3),
                },
            'dim': pb_def['dim'] - 1,
            'class': cb.CorrDim,
            'save_name': 'corrs_%s_pi%s' % (pb_def['name'], ch),
            'dump_variables': ['p' + ch],
            },
        })

    for ipm in ['p', 'm']:
        coefs.update({
            'E' + ipm + ch: {  # test+
                'requires': ['corrs_eta' + ch],
                'set_variables': [('corr_M', 'corrs_eta' + ch, 'pM')],
                'expression': 'ev_surface_integrate.iS.bYM%s(corr_M)' % ipm,
                'set_volume': 'bYp',
                'class': cb.CoefOne,
                },
            'F' + ipm + ch: {  # test+
                'requires': ['corrs_one' + ch, 'corrs_gamma_' + ipm],
                'set_variables': [('corr1_M', 'corrs_one' + ch, 'pM'),
                                  ('corr2_M', 'corrs_gamma_' + ipm, 'pM')],
                'expression': """dw_diffusion.iV.YM(mat2M.k, corr1_M, corr2_M)
                          - %e * ev_surface_integrate.iS.bYM%s(corr1_M)"""\
                                 % (1.0/param_h, ipm),
                'class': cb.CoefOne,
                },
            })

    for i_io in range(len(val['io_nd_grp'])):
        io = '%s_%d' % (ch, i_io + 1)

        coefs.update({
            'S' + io: { # [Rohan1] (4.28), test+
                'requires': ['corrs_gamma_' + io, 'pis_' + ch],
                'set_variables': set_corr_S,
                'expression': 'dw_diffusion.iV.Y%s(mat2%s.k,corr1_%s,corr2_%s)'\
                              % ((ch,) * 4),
                'dim': pb_def['dim'] - 1,
                'class': cb.CoefDim,
                },
            'P' + io: {  # test+
                'requires': ['pis_' + ch, 'corrs_pi' + ch],
                'set_variables': set_corr_cc,
                'expression': 'ev_surface_integrate.iS.bY%s(corr1_%s)'\
                              % (io, ch),
                'set_volume': 'bYp',
                'dim': pb_def['dim'] - 1,
                'class': cb.CoefDim,
                },
            'S_test' + io: {
                'requires': ['corrs_pi' + ch],
                'set_variables': [('corr1_' + ch, 'corrs_pi' + ch, 'p' + ch)],
                'expression':  '%e * ev_surface_integrate.iS.bY%s(corr1_%s)'\
                                % (1.0 / param_h, io, ch),
                'dim': pb_def['dim'] - 1,
                'class': cb.CoefDim,
                },
            })

        requirements.update({
            'corrs_gamma_' + io: {
                'requires': [],
                'variables': ['p' + ch, 'q' + ch],
                'ebcs': [],
                'epbcs': all_periodicY[ch],
                'lcbcs': ['imv' + ch],
                'equations': {
                    'eq_gamma': """dw_diffusion.iV.Y%s(mat2%s.k, q%s, p%s)
                                   + dw_volume_dot.iV.Y%s(q%s, ls%s)
                                   = %e * dw_surface_integrate.iS.bY%s(q%s)"""\
                                    % ((ch,) * 7 + (1.0/param_h, io, ch)),
                    'eq_imv': 'dw_volume_dot.iV.Y%s(lv%s, p%s) = 0'\
                              % ((ch,) * 3),
                    },
                'class': cb.CorrOne,
                'save_name': 'corrs_%s_gamma_%s' % (pb_def['name'], io),
                'dump_variables': ['p' + ch],
                },
            })

        for i_io2 in range(len(val['io_nd_grp'])):
            io2 = '%s_%d' % (ch, i_io2 + 1)
            io12 = '%s_%d' % (io, i_io2 + 1)
            coefs.update({
                'R' + io12: {  # test+
                    'requires': ['corrs_gamma_' + io2],
                    'set_variables': [('corr1_' + ch, 'corrs_gamma_' + io2,
                                       'p' + ch)],
                    'expression': 'ev_surface_integrate.iS.bY%s(corr1_%s)'\
                                  % (io, ch),
                    'set_volume': 'bYp',
                    'class': cb.CoefOne,
                    },
                })

###########################################################################
# solver, fe

solvers = {
    'ls': ('ls.scipy_direct', {}),
    'newton': ('nls.newton', {
        'i_max': 1,
    })
}

fe = {
    'chunk_size': 1000
}


# [Rohan1] Rohan E.: HOMOGENIZATION OF THE DARCY FLOW IN A DOUBLE-POROUS LAYER