Python Function.split方法代码示例

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
class TVPD(TV):
    """ Total variation using primal-dual Newton """

    def isPD(self): return True

    def updatePD(self):
        """ Update the parameters.
        parameters should be:
            - k(x) = factor inside TV
            - eps = regularization parameter
            - Vm = FunctionSpace for parameter. 
        ||f||_TV = int k(x) sqrt{|grad f|^2 + eps} dx
        """
        # primal dual variables
        self.Vw = FunctionSpace(self.Vm.mesh(), 'DG', 0)
        self.wH = Function(self.Vw*self.Vw)  # dual variable used in Hessian (re-scaled)
        #self.wH = nabla_grad(self.m)/sqrt(self.fTV) # full Hessian
        self.w = Function(self.Vw*self.Vw)  # dual variable for primal-dual, initialized at 0
        self.dm = Function(self.Vm)
        self.dw = Function(self.Vw*self.Vw)  
        self.testw = TestFunction(self.Vw*self.Vw)
        self.trialw = TrialFunction(self.Vw*self.Vw)
        # investigate convergence of dual variable
        self.dualres = self.w*sqrt(self.fTV) - nabla_grad(self.m)
        self.dualresnorm = inner(self.dualres, self.dualres)*dx
        self.normgraddm = inner(nabla_grad(self.dm), nabla_grad(self.dm))*dx
        # Hessian
        self.wkformPDhess = self.kovsq * ( \
        inner(nabla_grad(self.trial), nabla_grad(self.test)) - \
        0.5*( inner(self.wH, nabla_grad(self.test))*\
        inner(nabla_grad(self.trial), nabla_grad(self.m)) + \
        inner(nabla_grad(self.m), nabla_grad(self.test))*\
        inner(nabla_grad(self.trial), self.wH) ) / sqrt(self.fTV) \
        )*dx
        # update dual variable
        self.Mw = assemble(inner(self.trialw, self.testw)*dx)
        self.rhswwk = inner(-self.w, self.testw)*dx + \
        inner(nabla_grad(self.m)+nabla_grad(self.dm), self.testw) \
        /sqrt(self.fTV)*dx + \
        inner(-inner(nabla_grad(self.m),nabla_grad(self.dm))* \
        self.wH/sqrt(self.fTV), self.testw)*dx

    def compute_dw(self, dm):
        """ Compute dw """
        setfct(self.dm, dm)
        b = assemble(self.rhswwk)
        solve(self.Mw, self.dw.vector(), b)


    def update_w(self, alpha):
        """ update w and re-scale wH """
        self.w.vector().axpy(alpha, self.dw.vector())
        # project each w (coord-wise) onto unit sphere to get wH
        (wx, wy) = self.w.split(deepcopy=True)
        wxa, wya = wx.vector().array(), wy.vector().array()
        normw = np.sqrt(wxa**2 + wya**2)
        factorw = [max(1.0, ii) for ii in normw]
        setfct(wx, wxa/factorw)
        setfct(wy, wya/factorw)
        assign(self.wH.sub(0), wx)
        assign(self.wH.sub(1), wy)
        # check
        (wx,wy) = self.wH.split(deepcopy=True)
        wxa, wya = wx.vector().array(), wy.vector().array()
        assert np.amax(np.sqrt(wxa**2 + wya**2)) <= 1.0 + 1e-14
        # Print results
        dualresnorm = assemble(self.dualresnorm)
        normgraddm = assemble(self.normgraddm)
        print 'line search dual variable: max(|w|)={}, err(w,df)={}, |grad(dm)|={}'.\
        format(np.amax(np.sqrt(normw)), np.sqrt(dualresnorm), np.sqrt(normgraddm))

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
    def solve(self, problem):
        self.problem = problem
        doSave = problem.doSave
        save_this_step = False
        onlyVel = problem.saveOnlyVel
        dt = self.metadata['dt']

        nu = Constant(self.problem.nu)
        self.tc.init_watch('init', 'Initialization', True, count_to_percent=False)
        self.tc.init_watch('solve', 'Running nonlinear solver', True, count_to_percent=True)
        self.tc.init_watch('next', 'Next step assignments', True, count_to_percent=True)
        self.tc.init_watch('saveVel', 'Saved velocity', True)

        self.tc.start('init')

        mesh = self.problem.mesh

        # Define function spaces (P2-P1)
        self.V = VectorFunctionSpace(mesh, "Lagrange", 2)  # velocity
        self.Q = FunctionSpace(mesh, "Lagrange", 1)  # pressure
        self.W = MixedFunctionSpace([self.V, self.Q])
        self.PS = FunctionSpace(mesh, "Lagrange", 2)  # partial solution (must be same order as V)
        self.D = FunctionSpace(mesh, "Lagrange", 1)   # velocity divergence space

        # to assign solution in space W.sub(0) to Function(V) we need FunctionAssigner (cannot be assigned directly)
        fa = FunctionAssigner(self.V, self.W.sub(0))
        velSp = Function(self.V)

        problem.initialize(self.V, self.Q, self.PS, self.D)

        # Define unknown and test function(s) NS
        v, q = TestFunctions(self.W)
        w = Function(self.W)
        dw = TrialFunction(self.W)
        u, p = split(w)

        # Define fields
        n = FacetNormal(mesh)
        I = Identity(u.geometric_dimension())
        theta = 0.5  # Crank-Nicholson
        k = Constant(self.metadata['dt'])

        # Initial conditions: u0 velocity at previous time step u1 velocity two time steps back p0 previous pressure
        [u0, p0] = self.problem.get_initial_conditions([{'type': 'v', 'time': 0.0}, {'type': 'p', 'time': 0.0}])

        if doSave:
            problem.save_vel(False, u0, 0.0)

        # boundary conditions
        bcu, bcp = problem.get_boundary_conditions(self.bc == 'outflow', self.W.sub(0), self.W.sub(1))
        # NT bcp is not used

        # Define steady part of the equation
        def T(u):
            return -p * I + 2.0 * nu * sym(grad(u))

        def F(u, v, q):
            return (inner(T(u), grad(v)) - q * div(u)) * dx + inner(grad(u) * u, v) * dx

        # Define variational forms
        F_ns = (inner((u - u0), v) / k) * dx + (1.0 - theta) * F(u0, v, q) + theta * F(u, v, q)
        J_ns = derivative(F_ns, w, dw)
        # J_ns = derivative(F_ns, w)  # did not work

        # NS_problem = NonlinearVariationalProblem(F_ns, w, bcu, J_ns, form_compiler_parameters=ffc_options)
        NS_problem = NonlinearVariationalProblem(F_ns, w, bcu, J_ns)
        # (var. formulation, unknown, Dir. BC, jacobian, optional)
        NS_solver = NonlinearVariationalSolver(NS_problem)

        prm = NS_solver.parameters
        prm['newton_solver']['absolute_tolerance'] = 1E-08
        prm['newton_solver']['relative_tolerance'] = 1E-08
        # prm['newton_solver']['maximum_iterations'] = 45
        # prm['newton_solver']['relaxation_parameter'] = 1.0
        prm['newton_solver']['linear_solver'] = 'mumps'

        info(NS_solver.parameters, True)

        self.tc.end('init')

        # Time-stepping
        info("Running of direct method")
        ttime = self.metadata['time']
        t = dt
        step = 1
        while t < (ttime + dt/2.0):
            info("t = %f" % t)
            self.problem.update_time(t, step)
            if self.MPI_rank == 0:
                problem.write_status_file(t)

            if doSave:
                save_this_step = problem.save_this_step

            # Compute
            begin("Solving NS ....")
            try:
                self.tc.start('solve')
                NS_solver.solve()
                self.tc.end('solve')
#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
    def _evaluateGlobalMixedEstimator(cls, mu, w, coeff_field, pde, f, quadrature_degree, vectorspace_type='BDM'):
        """Evaluation of global mixed equilibrated estimator."""
        # set quadrature degree
#        quadrature_degree_old = parameters["form_compiler"]["quadrature_degree"]
#        parameters["form_compiler"]["quadrature_degree"] = quadrature_degree
#        logger.debug("residual quadrature order = " + str(quadrature_degree))

        # prepare numerical flux and f
        sigma_mu, f_mu = evaluate_numerical_flux(w, mu, coeff_field, f)

        # ###################
        # ## MIXED PROBLEM ##
        # ###################

        # get setup data for mixed problem
        V = w[mu]._fefunc.function_space()
        mesh = V.mesh()
        degree = element_degree(w[mu]._fefunc)

        # create function spaces
        DG0 = FunctionSpace(mesh, 'DG', 0)
        DG0_dofs = [DG0.dofmap().cell_dofs(c.index())[0] for c in cells(mesh)]
        RT = FunctionSpace(mesh, vectorspace_type, degree)
        W = RT * DG0

        # setup boundary conditions
#        bcs = pde.create_dirichlet_bcs(W.sub(1))

        # debug ===
        # from dolfin import DOLFIN_EPS, DirichletBC
        # def boundary(x):
        #     return x[0] < DOLFIN_EPS or x[0] > 1.0 + DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 + DOLFIN_EPS
        # bcs = [DirichletBC(W.sub(1), Constant(0.0), boundary)]
        # === debug

        # create trial and test functions
        (sigma, u) = TrialFunctions(W)
        (tau, v) = TestFunctions(W)

        # define variational form
        a_eq = (dot(sigma, tau) + div(tau) * u + div(sigma) * v) * dx
        L_eq = (- f_mu * v + dot(sigma_mu, tau)) * dx

        # compute solution
        w_eq = Function(W)
        solve(a_eq == L_eq, w_eq)
        (sigma_mixed, u_mixed) = w_eq.split()

        # #############################
        # ## EQUILIBRATION ESTIMATOR ##
        # #############################

        # evaluate error estimator
        dg0 = TestFunction(DG0)
        eta_mu = inner(sigma_mu, sigma_mu) * dg0 * dx
        eta_T = assemble(eta_mu, form_compiler_parameters={'quadrature_degree': quadrature_degree})
        eta_T = np.array([sqrt(e) for e in eta_T])

        # evaluate global error
        eta = sqrt(sum(i**2 for i in eta_T))
        # reorder array entries for local estimators
        eta_T = eta_T[DG0_dofs]

        # restore quadrature degree
#        parameters["form_compiler"]["quadrature_degree"] = quadrature_degree_old

        return eta, FlatVector(eta_T)

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
def solve_fixed_point(
        mesh,
        W_element, P_element, Q_element,
        u0, p0, theta0,
        kappa, rho, mu, cp,
        g, extra_force,
        heat_source,
        u_bcs, p_bcs,
        theta_dirichlet_bcs,
        theta_neumann_bcs,
        my_dx, my_ds,
        max_iter,
        tol
        ):
    # Solve the coupled heat-Stokes equation approximately. Do this
    # iteratively by solving the heat equation, then solving Stokes with the
    # updated heat, the heat equation with the updated velocity and so forth
    # until the change is 'small'.
    WP = FunctionSpace(mesh, MixedElement([W_element, P_element]))
    Q = FunctionSpace(mesh, Q_element)
    # Initialize functions.
    up0 = Function(WP)
    u0, p0 = up0.split()

    theta1 = Function(Q)
    for _ in range(max_iter):
        heat_problem = heat.Heat(
            Q,
            kappa=kappa,
            rho=rho(theta0),
            cp=cp,
            convection=u0,
            source=heat_source,
            dirichlet_bcs=theta_dirichlet_bcs,
            neumann_bcs=theta_neumann_bcs,
            my_dx=my_dx,
            my_ds=my_ds
            )

        theta1.assign(heat_problem.solve_stationary())

        # Solve problem for velocity, pressure.
        f = rho(theta0) * g  # coupling
        if extra_force:
            f += as_vector((extra_force[0], extra_force[1], 0.0))
        # up1 = up0.copy()
        stokes.stokes_solve(
            up0,
            mu,
            u_bcs, p_bcs,
            f,
            my_dx=my_dx,
            tol=1.0e-10,
            verbose=False,
            maxiter=1000
            )

        # from dolfin import plot
        # plot(u0)
        # plot(theta0)

        theta_diff = errornorm(theta0, theta1)
        info('||theta - theta0|| = {:e}'.format(theta_diff))
        # info('||u - u0||         = {:e}'.format(u_diff))
        # info('||p - p0||         = {:e}'.format(p_diff))
        # diff = theta_diff + u_diff + p_diff
        diff = theta_diff
        info('sum = {:e}'.format(diff))

        # # Show the iterates.
        # plot(theta0, title='theta0')
        # plot(u0, title='u0')
        # interactive()
        # #exit()
        if diff < tol:
            break

        theta0.assign(theta1)

    # Create a *deep* copy of u0, p0, to be able to deal with them as
    # actually separate entities.
    u0, p0 = up0.split(deepcopy=True)
    return u0, p0, theta0

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
    def ab2tr_step0(u0,
                    P,
                    f,  # right-hand side
                    rho,
                    mu,
                    dudt_bcs=[],
                    p_bcs=[],
                    eps=1.0e-4,  # relative error tolerance
                    verbose=True
                    ):
        # Make sure that the initial velocity is divergence-free.
        alpha = norm(u0, 'Hdiv0')
        if abs(alpha) > DOLFIN_EPS:
            warn('Initial velocity not divergence-free (||u||_div = %e).'
                 % alpha
                 )
        # Get the initial u0' and p0 by solving the linear equation system
        #
        #     [M   C] [u0']   [f0 - (K+N(u0)u0)]
        #     [C^T 0] [p0 ] = [ g0'            ],
        #
        # i.e.,
        #
        #     rho u0' + nabla(p0) = f0 + mu\Delta(u0) - rho u0.nabla(u0),
        #     div(u0')            = 0.
        #
        W = u0.function_space()
        WP = W*P

        # Translate the boundary conditions into product space. See
        # <http://fenicsproject.org/qa/703/boundary-conditions-in-product-space>.
        dudt_bcs_new = []
        for dudt_bc in dudt_bcs:
            dudt_bcs_new.append(DirichletBC(WP.sub(0),
                                            dudt_bc.value(),
                                            dudt_bc.user_sub_domain()))
        p_bcs_new = []
        for p_bc in p_bcs:
            p_bcs_new.append(DirichletBC(WP.sub(1),
                                         p_bc.value(),
                                         p_bc.user_sub_domain()))

        new_bcs = dudt_bcs_new + p_bcs_new

        (u, p) = TrialFunctions(WP)
        (v, q) = TestFunctions(WP)

        #a = rho * dot(u, v) * dx + dot(grad(p), v) * dx \
        a = rho * inner(u, v) * dx - p * div(v) * dx \
            - div(u) * q * dx
        L = _rhs_weak(u0, v, f, rho, mu)

        A, b = assemble_system(a, L, new_bcs)

        # Similar preconditioner as for the Stokes problem.
        # TODO implement something better!
        prec = rho * inner(u, v) * dx \
            - p*q*dx
        M, _ = assemble_system(prec, L, new_bcs)

        solver = KrylovSolver('gmres', 'amg')

        solver.parameters['monitor_convergence'] = verbose
        solver.parameters['report'] = verbose
        solver.parameters['absolute_tolerance'] = 0.0
        solver.parameters['relative_tolerance'] = 1.0e-6
        solver.parameters['maximum_iterations'] = 10000

        # Associate operator (A) and preconditioner matrix (M)
        solver.set_operators(A, M)
        #solver.set_operator(A)

        # Solve
        up = Function(WP)
        solver.solve(up.vector(), b)

        # Get sub-functions
        dudt0, p0 = up.split()

        # Choosing the first step size for the trapezoidal rule can be tricky.
        # Chapters 2.7.4a, 2.7.4e of the book
        #
        #     Incompressible flow and the finite element method,
        #     volume 1: advection-diffusion;
        #     P.M. Gresho, R.L. Sani,
        #
        # give some hints.
        #
        #     eps ... relative error tolerance
        #     tau ... estimate of the initial 'time constant'
        tau = None
        if tau:
            dt0 = tau * eps**(1.0/3.0)
        else:
            # Choose something 'reasonably small'.
            dt0 = 1.0e-3
        # Alternative:
        # Use a dissipative scheme like backward Euler or BDF2 for the first
        # couple of steps. This makes sure that noisy initial data is damped
        # out.
#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
def solve(
        mesh,
        W_element, P_element, Q_element,
        u0, p0, theta0,
        kappa, rho, mu, cp,
        g, extra_force,
        heat_source,
        u_bcs, p_bcs,
        theta_dirichlet_bcs,
        theta_neumann_bcs,
        dx_submesh, ds_submesh
        ):
    # First do a fixed_point iteration. This is usually quite robust and leads
    # to a point from where Newton can converge reliably.
    u0, p0, theta0 = solve_fixed_point(
        mesh,
        W_element, P_element, Q_element,
        theta0,
        kappa, rho, mu, cp,
        g, extra_force,
        heat_source,
        u_bcs, p_bcs,
        theta_dirichlet_bcs,
        theta_neumann_bcs,
        my_dx=dx_submesh,
        my_ds=ds_submesh,
        max_iter=100,
        tol=1.0e-8
        )

    WPQ = FunctionSpace(
        mesh, MixedElement([W_element, P_element, Q_element])
        )
    uptheta0 = Function(WPQ)

    # Initial guess
    assign(uptheta0.sub(0), u0)
    assign(uptheta0.sub(1), p0)
    assign(uptheta0.sub(2), theta0)

    rho_const = rho(_average(theta0))

    stokes_heat_problem = StokesHeat(
        WPQ,
        kappa, rho, rho_const, mu, cp,
        g, extra_force,
        heat_source,
        u_bcs, p_bcs,
        theta_dirichlet_bcs=theta_dirichlet_bcs,
        theta_neumann_bcs=theta_neumann_bcs,
        my_dx=dx_submesh,
        my_ds=ds_submesh
        )

    # solver = FixedPointSolver()
    from dolfin import PETScSNESSolver
    solver = PETScSNESSolver()
    # http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESType.html
    solver.parameters['method'] = 'newtonls'
    # The Jacobian system for Stokes (+heat) are hard to solve.
    # Use LU for now.
    solver.parameters['linear_solver'] = 'lu'
    solver.parameters['maximum_iterations'] = 100
    # TODO tighten tolerance. might not always work though...
    solver.parameters['absolute_tolerance'] = 1.0e-3
    solver.parameters['relative_tolerance'] = 0.0
    solver.parameters['report'] = True

    solver.solve(stokes_heat_problem, uptheta0.vector())

    # u0, p0, theta0 = split(uptheta0)
    # Create a *deep* copy of u0, p0, theta0 to be able to deal with them
    # as actually separate entities vectors.
    u0, p0, theta0 = uptheta0.split(deepcopy=True)
    return u0, p0, theta0

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
def solve(W, P,
          mu,
          u_bcs, p_bcs,
          f,
          verbose=True,
          tol=1.0e-10
          ):
    # Some initial sanity checks.
    assert mu > 0.0

    WP = MixedFunctionSpace([W, P])

    # Translate the boundary conditions into the product space.
    # This conditional loop is able to deal with conditions of the kind
    #
    #     DirichletBC(W.sub(1), 0.0, right_boundary)
    #
    new_bcs = []
    for k, bcs in enumerate([u_bcs, p_bcs]):
        for bc in bcs:
            space = bc.function_space()
            C = space.component()
            if len(C) == 0:
                new_bcs.append(DirichletBC(WP.sub(k),
                                           bc.value(),
                                           bc.domain_args[0]))
            elif len(C) == 1:
                new_bcs.append(DirichletBC(WP.sub(k).sub(int(C[0])),
                                           bc.value(),
                                           bc.domain_args[0]))
            else:
                raise RuntimeError('Illegal number of subspace components.')

    # Define variational problem
    (u, p) = TrialFunctions(WP)
    (v, q) = TestFunctions(WP)

    # Build system.
    # The sign of the div(u)-term is somewhat arbitrary since the right-hand
    # side is 0 here. We can either make the system symmetric or positive-
    # definite.
    # On a second note, we have
    #
    #    \int grad(p).v = - \int p * div(v) + \int_\Gamma p n.v.
    #
    # Since, we have either p=0 or n.v=0 on the boundary, we could as well
    # replace the term dot(grad(p), v) by -p*div(v).
    #
    a = mu * inner(grad(u), grad(v))*dx \
      - p * div(v) * dx \
      - q * div(u) * dx
    #a = mu * inner(grad(u), grad(v))*dx + dot(grad(p), v) * dx \
    #  - div(u) * q * dx
    L = dot(f, v)*dx
    A, b = assemble_system(a, L, new_bcs)

    if has_petsc():
        # For an assortment of preconditioners, see
        #
        #     Performance and analysis of saddle point preconditioners
        #     for the discrete steady-state Navier-Stokes equations;
        #     H.C. Elman, D.J. Silvester, A.J. Wathen;
        #     Numer. Math. (2002) 90: 665-688;
        #     <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.145.3554>.
        #
        # Set up field split.
        W = SubSpace(WP, 0)
        P = SubSpace(WP, 1)
        u_dofs = W.dofmap().dofs()
        p_dofs = P.dofmap().dofs()
        prec = PETScPreconditioner()
        prec.set_fieldsplit([u_dofs, p_dofs], ['u', 'p'])

        PETScOptions.set('pc_type', 'fieldsplit')
        PETScOptions.set('pc_fieldsplit_type', 'additive')
        PETScOptions.set('fieldsplit_u_pc_type', 'lu')
        PETScOptions.set('fieldsplit_p_pc_type', 'jacobi')

        ## <http://scicomp.stackexchange.com/questions/7288/which-preconditioners-and-solver-in-petsc-for-indefinite-symmetric-systems-sho>
        #PETScOptions.set('pc_type', 'fieldsplit')
        ##PETScOptions.set('pc_fieldsplit_type', 'schur')
        ##PETScOptions.set('pc_fieldsplit_schur_fact_type', 'upper')
        #PETScOptions.set('pc_fieldsplit_detect_saddle_point')
        ##PETScOptions.set('fieldsplit_u_pc_type', 'lsc')
        ##PETScOptions.set('fieldsplit_u_ksp_type', 'preonly')

        #PETScOptions.set('pc_type', 'fieldsplit')
        #PETScOptions.set('fieldsplit_u_pc_type', 'hypre')
        #PETScOptions.set('fieldsplit_u_ksp_type', 'preonly')
        #PETScOptions.set('fieldsplit_p_pc_type', 'jacobi')
        #PETScOptions.set('fieldsplit_p_ksp_type', 'preonly')

        ## From PETSc/src/ksp/ksp/examples/tutorials/ex42-fsschur.opts:
        #PETScOptions.set('pc_type', 'fieldsplit')
        #PETScOptions.set('pc_fieldsplit_type', 'SCHUR')
        #PETScOptions.set('pc_fieldsplit_schur_fact_type', 'UPPER')
        #PETScOptions.set('fieldsplit_p_ksp_type', 'preonly')
        #PETScOptions.set('fieldsplit_u_pc_type', 'bjacobi')

        ## From
#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
def RunJob(Tb, mu_value, path):
    runtimeInit = clock()

    tfile = File(path + '/t6t.pvd')
    mufile = File(path + "/mu.pvd")
    ufile = File(path + '/velocity.pvd')
    gradpfile = File(path + '/gradp.pvd')
    pfile = File(path + '/pstar.pvd')
    parameters = open(path + '/parameters', 'w', 0)
    vmeltfile = File(path + '/vmelt.pvd')
    rhofile = File(path + '/rhosolid.pvd')

    for name in dir():
        ev = str(eval(name))
        if name[0] != '_' and ev[0] != '<':
            parameters.write(name + ' = ' + ev + '\n')

    temp_values = [27. + 273, Tb + 273, 1300. + 273, 1305. + 273]
    dTemp = temp_values[3] - temp_values[0]
    temp_values = [x / dTemp for x in temp_values]  # non dimensionalising temp

    mu_a = mu_value  # this was taken from the blankenbach paper, can change..
    
    Ep = b / dTemp

    mu_bot = exp(-Ep * (temp_values[3] * dTemp - 1573) + cc) * mu_a

    Ra = rho_0 * alpha * g * dTemp * h**3 / (kappa_0 * mu_a)
    w0 = rho_0 * alpha * g * dTemp * h**2 / mu_a
    tau = h / w0
    p0 = mu_a * w0 / h

    print(mu_a, mu_bot, Ra, w0, p0)

    vslipx = 1.6e-09 / w0
    vslip = Constant((vslipx, 0.0))  # nondimensional
    noslip = Constant((0.0, 0.0))

    dt = 3.E11 / tau
    tEnd = 3.E13 / tau  # non-dimensionalising times

    class PeriodicBoundary(SubDomain):
        def inside(self, x, on_boundary):
            return left(x, on_boundary)

        def map(self, x, y):
            y[0] = x[0] - MeshWidth
            y[1] = x[1]

    pbc = PeriodicBoundary()

    class TempExp(Expression):
        def eval(self, value, x):
            if x[1] >= LAB(x):
                value[0] = temp_values[0] + (temp_values[1] - temp_values[0]) * (MeshHeight - x[1]) / (MeshHeight - LAB(x))
            else:
                value[0] = temp_values[3] - (temp_values[3] - temp_values[2]) * (x[1]) / (LAB(x))

    class FluidTemp(Expression):
        def eval(self, value, x):
            if value[0] < 1295:
                value[0] = 1295

    mesh = RectangleMesh(Point(0.0, 0.0), Point(MeshWidth, MeshHeight), nx, ny)

    Svel = VectorFunctionSpace(mesh, 'CG', 2, constrained_domain=pbc)
    Spre = FunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)
    Stemp = FunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)
    Smu = FunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)
    Sgradp = VectorFunctionSpace(mesh, 'CG', 2, constrained_domain=pbc)
    Srho = FunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)
    S0 = MixedFunctionSpace([Svel, Spre, Stemp])

    u = Function(S0)
    v, p, T = split(u)
    v_t, p_t, T_t = TestFunctions(S0)

    T0 = interpolate(TempExp(), Stemp)

    muExp = Expression('exp(-Ep * (T_val * dTemp - 1573) + cc * x[2] / meshHeight)', Smu.ufl_element(),
                        Ep=Ep, dTemp=dTemp, cc=cc, meshHeight=MeshHeight, T_val=T0)

    mu = interpolate(muExp, Smu)

    rhosolid = Function(Srho)
    deltarho = Function(Srho)

    v0 = Function(Svel)
    vmelt = Function(Svel)

    v_theta = (1. - theta)*v0 + theta*v

    T_theta = (1. - theta)*T + theta*T0

    r_v = (inner(sym(grad(v_t)), 2.*mu*sym(grad(v))) \
        - div(v_t)*p \
        - T*v_t[1] )*dx

    r_p = p_t*div(v)*dx

#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
class NavierStokesCylindrical(NonlinearProblem):

    def __init__(self, WP, rho, mu, f, u_bcs, p_bcs, dx, stabilization=None):
        super(NavierStokesCylindrical, self).__init__()
        #
        self.WP = WP
        # Translate the boundary conditions into the product space.
        self.bcs = []
        for k, bcs in enumerate([u_bcs, p_bcs]):
            for bc in bcs:
                space = bc.function_space()
                C = space.component()
                if len(C) == 0:
                    self.bcs.append(DirichletBC(self.WP.sub(k),
                                                bc.value(),
                                                bc.domain_args[0]))
                elif len(C) == 1:
                    self.bcs.append(DirichletBC(self.WP.sub(k).sub(int(C[0])),
                                                bc.value(),
                                                bc.domain_args[0]))
                else:
                    raise RuntimeError('Illegal number of '
                                       'subspace components.'
                                       )
        #
        self.up = Function(self.WP)
        u, p = self.up.split()
        v, q = TestFunctions(self.WP)
        # Momentum equation.
        self.F0 = momentum_equation(u, v,
                                    p,
                                    f,
                                    rho, mu,
                                    stabilization=stabilization,
                                    dx=dx
                                    )
        # div_u = 1/r * div(r*u)
        r = Expression('x[0]', degree=1, domain=WP.mesh())
        self.F0 += (1.0 / r * (r * u[0]).dx(0) + u[1].dx(1)) * q \
            * 2 * pi * r * dx

        self.jacobian = derivative(self.F0, self.up)
        self.reset_sparsity = True
        return

    def F(self, b, x):
        self.up.vector()[:] = x
        assemble(self.F0,
                 tensor=b,
                 form_compiler_parameters={'optimize': True}
                 )
        for bc in self.bcs:
            bc.apply(b, x)
        return

    def J(self, A, x):
        self.up.vector()[:] = x
        assemble(self.jacobian,
                 tensor=A,
                 reset_sparsity=self.reset_sparsity,
                 form_compiler_parameters={'optimize': True}
                 )
        for bc in self.bcs:
            bc.apply(A)
        self.reset_sparsity = False
        return

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import split [as 别名]
def run_with_params(Tb, mu_value, k_s, path):
    run_time_init = clock()

    mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(mesh_width, mesh_width, mesh_height), nx, ny, nz)

    pbc = PeriodicBoundary()

    WE = VectorElement('CG', mesh.ufl_cell(), 2)
    SE = FiniteElement('CG', mesh.ufl_cell(), 1)
    WSSS = FunctionSpace(mesh, MixedElement(WE, SE, SE, SE), constrained_domain=pbc)
    # W = FunctionSpace(mesh, WE, constrained_domain=pbc)
    # S = FunctionSpace(mesh, SE, constrained_domain=pbc)
    W = WSSS.sub(0).collapse()
    S = WSSS.sub(1).collapse()

    temperature_vals = [27.0 + 273, Tb + 273, 1300.0 + 273, 1305.0 + 273]
    temp_prof = TemperatureProfile(temperature_vals, element=S.ufl_element())

    mu_a = mu_value  # this was taken from the Blankenbach paper, can change

    Ep = b / temp_prof.delta

    mu_bot = exp(-Ep * (temp_prof.bottom * temp_prof.delta - 1573.0) + cc) * mu_a

    # TODO: verify exponentiation
    Ra = rho_0 * alpha * g * temp_prof.delta * h ** 3 / (kappa_0 * mu_a)
    w0 = rho_0 * alpha * g * temp_prof.delta * h ** 2 / mu_a
    tau = h / w0
    p0 = mu_a * w0 / h

    log(mu_a, mu_bot, Ra, w0, p0)

    slip_vx = 1.6E-09 / w0  # Non-dimensional
    slip_velocity = Constant((slip_vx, 0.0, 0.0))
    zero_slip = Constant((0.0, 0.0, 0.0))

    time_step = 3.0E11 / tau * 2

    dt = Constant(time_step)
    t_end = 3.0E15 / tau / 5.0  # Non-dimensional times

    u = Function(WSSS)

    # Instead of TrialFunctions, we use split(u) for our non-linear problem
    v, p, T, Tf = split(u)
    v_t, p_t, T_t, Tf_t = TestFunctions(WSSS)

    T0 = interpolate(temp_prof, S)

    mu_exp = Expression('exp(-Ep * (T_val * dTemp - 1573.0) + cc * x[2] / mesh_height)',
                       Ep=Ep, dTemp=temp_prof.delta, cc=cc, mesh_height=mesh_height, T_val=T0,
                       element=S.ufl_element())

    Tf0 = interpolate(temp_prof, S)

    mu = Function(S)
    v0 = Function(W)

    v_theta = (1.0 - theta) * v0 + theta * v

    T_theta = (1.0 - theta) * T0 + theta * T

    Tf_theta = (1.0 - theta) * Tf0 + theta * Tf

    # TODO: Verify forms

    r_v = (inner(sym(grad(v_t)), 2.0 * mu * sym(grad(v)))
           - div(v_t) * p
           - T * v_t[2]) * dx

    r_p = p_t * div(v) * dx

    heat_transfer = Constant(k_s) * (Tf_theta - T_theta) * dt

    r_T = (T_t * ((T - T0) + dt * inner(v_theta, grad(T_theta)))  # TODO: Inner vs dot
           + (dt / Ra) * inner(grad(T_t), grad(T_theta))
           - T_t * heat_transfer) * dx

    v_melt = Function(W)
    z_hat = Constant((0.0, 0.0, 1.0))

    # TODO: inner -> dot, take out Tf_t
    r_Tf = (Tf_t * ((Tf - Tf0) + dt * inner(v_melt, grad(Tf_theta)))
            + Tf_t * heat_transfer) * dx

    r = r_v + r_p + r_T + r_Tf

    bcv0 = DirichletBC(WSSS.sub(0), zero_slip, top)
    bcv1 = DirichletBC(WSSS.sub(0), slip_velocity, bottom)
    bcv2 = DirichletBC(WSSS.sub(0).sub(1), Constant(0.0), back)
    bcv3 = DirichletBC(WSSS.sub(0).sub(1), Constant(0.0), front)

    bcp0 = DirichletBC(WSSS.sub(1), Constant(0.0), bottom)
    bct0 = DirichletBC(WSSS.sub(2), Constant(temp_prof.surface), top)
    bct1 = DirichletBC(WSSS.sub(2), Constant(temp_prof.bottom), bottom)
    bctf1 = DirichletBC(WSSS.sub(3), Constant(temp_prof.bottom), bottom)

    bcs = [bcv0, bcv1, bcv2, bcv3, bcp0, bct0, bct1, bctf1]

    t = 0
#.........这里部分代码省略.........