I attach the codes here:

from __future__ import absolute_import from sfepy.examples.dg.example_dg_common import * from sfepy import data_dir import numpy as nm from os.path import join as pjoin from sfepy.discrete.fem import Mesh from sfepy.discrete.fem.meshio import UserMeshIO import matplotlib.pyplot as plt from sfepy.discrete import (FieldVariable, Material, Integral, Function, Equation, Equations, Problem) from sfepy.discrete.problem import prepare_matrix from sfepy.discrete.fem import Mesh, FEDomain, Field from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC, InitialCondition from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton from sfepy.solvers.ts_solvers import SimpleTimeSteppingSolver from sfepy.base.base import IndexedStruct from sfepy.base.base import Struct, configure_output, output

def main(): outputs_folder = "../outputs"

problem_name = "Fluid_dynamics"
output_folder = pjoin(outputs_folder, problem_name)
output_format = "vtk"
save_timestn = 100
clear_folder(pjoin(output_folder, "*." + output_format))
configure_output({'output_screen': True,
                  'output_log_name':
                   pjoin(output_folder,
                         f"last_run_{problem_name}.txt")})
						 
# Geometry meshing.
filename_mesh=None
X1 = 0
XN = 1.
n_nod = 1001
n_el = n_nod - 1

if filename_mesh is None:
    filename_mesh = get_gen_1D_mesh_hook(X1, XN, n_nod)

# Domain, variable setting.
domain = FEDomain('domain', filename_mesh)
omega = domain.create_region('Omega', 'all')

left = domain.create_region('Gamma_Left',
                            'vertices in x < 0.00001', 'facet')
right = domain.create_region('Gamma_Right',
                             'vertices in x > 0.99999', 'facet')

field = Field.from_args('f', nm.float64, 'scalar', omega,
                        approx_order=2)
rhow = FieldVariable('rhow', 'unknown', field, history=1)
v = FieldVariable('v', 'test', field, primary_var_name='rhow')

porosity = Material('porosity', val=0.3)

def get_diffusivity(ts, coors, problem, equations=None, mode=None, **kwargs):
    if mode == 'qp':
    # rho-field values in quadrature points coordinates given by integral i
    # - they are the same as in `coors` argument.
        rhow_values = problem.evaluate('ev_volume_integrate.i.Omega(rhow)',
                                mode='qp', verbose=False)
        val = 0.1*rhow_values
        output('diffusivity: min:', val.min(), 'max:', val.max())
        val.shape = (coors.shape[0], 1, 1)
return {'val' : val}

diffusivity = Material('diffusivity' , function=Function('get_diffusivity', get_diffusivity))

integral = Integral('i', order=2*2)

MassT = Term.new('dw_volume_dot.i.Omega(porosity.val, v, drhow/dt)', integral, omega, porosity=porosity, v=v, rhow=rhow)
StiffT = Term.new('dw_laplace.i.Omega(diffusivity.val, v, rhow)', integral, omega, diffusivity=diffusivity, v=v, rhow=rhow)

eq = Equation('balance', MassT + StiffT)
eqs = Equations([eq])

# Boundary conditions.
ebc1 = EssentialBC('Dens1', left, {'rhow.0' : 0.0})
ebc2 = EssentialBC('Dens2', right, {'rhow.0' : 1.0})

# Initial conditions.
def get_ic(coors, ic):
    return coors[:,0]
ic_fun = Function('ic_fun', get_ic)
ic = InitialCondition('ic', omega, {'rhow.0' : ic_fun})

# Create Problem.
pb = Problem(problem_name,
             equations=eqs,
             conf=Struct(options={"save_times": save_timestn}))
pb.set_bcs(ebcs=Conditions([ebc1, ebc2]))
pb.setup_output(output_dir=output_folder, output_format=output_format)
pb.set_ics(Conditions([ic]))
pb.update_materials()


ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({}, lin_solver=ls, status=nls_status)
tss = SimpleTimeSteppingSolver({'t0' : 0.0, 't1' : 1.0, 'n_step' : 100},
                               nls=nls, context=pb, verbose=True)
							   
# ---------
# | Solve |
# ---------
pb.set_solver(tss)
state_end = pb.solve()

# ----------
# | Plot 1D|
# ----------
load_and_plot_fun(output_folder, "domain",
                  t0, t1, min(tn, save_timestn),
                  ic_fun)

if __name__ == '__main__': main()

thanks, it works, though still existing some small bugs