Python Function.dx方法代码示例

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import dx [as 别名]
def compute_velocity_correction(
    ui, p0, p1, u_bcs, rho, mu, dt, rotational_form, my_dx, tol, verbose
):
    """Compute the velocity correction according to

    .. math::

        U = u_0 - \\frac{dt}{\\rho} \\nabla (p_1-p_0).
    """
    W = ui.function_space()
    P = p1.function_space()

    u = TrialFunction(W)
    v = TestFunction(W)
    a3 = dot(u, v) * my_dx
    phi = Function(P)
    phi.assign(p1)
    if p0:
        phi -= p0
    if rotational_form:
        r = SpatialCoordinate(W.mesh())[0]
        div_ui = 1 / r * (r * ui[0]).dx(0) + ui[1].dx(1)
        phi += mu * div_ui
    L3 = dot(ui, v) * my_dx - dt / rho * (phi.dx(0) * v[0] + phi.dx(1) * v[1]) * my_dx
    u1 = Function(W)
    solve(
        a3 == L3,
        u1,
        bcs=u_bcs,
        solver_parameters={
            "linear_solver": "iterative",
            "symmetric": True,
            "preconditioner": "hypre_amg",
            "krylov_solver": {
                "relative_tolerance": tol,
                "absolute_tolerance": 0.0,
                "maximum_iterations": 100,
                "monitor_convergence": verbose,
            },
        },
    )
    # u = project(ui - k/rho * grad(phi), V)
    # div_u = 1/r * div(r*u)
    r = SpatialCoordinate(W.mesh())[0]
    div_u1 = 1.0 / r * (r * u1[0]).dx(0) + u1[1].dx(1)
    info("||u||_div = {:e}".format(sqrt(assemble(div_u1 * div_u1 * my_dx))))
    return u1

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import dx [as 别名]
ds = Measure("ds", subdomain_data=boundary_parts)

# Define forms (dont redefine functions used here)
# step 1
u0 = Function(V)
u1 = Function(V)
p0 = Function(Q)
u_tent = TrialFunction(V)
v = TestFunction(V)
# U_ = 1.5*u0 - 0.5*u1
# nonlinearity = inner(dot(0.5 * (u_tent.dx(0) + u0.dx(0)), U_), v) * dx
# F_tent = (1./dt)*inner(u_tent - u0, v) * dx + nonlinearity\
#     + nu*inner(0.5 * (u_tent.dx(0) + u0.dx(0)), v.dx(0)) * dx + inner(p0.dx(0), v) * dx\
#     - inner(f, v)*dx     # solve to u_
# using explicite scheme: so LHS has interpretation as heat equation, RHS are sources
F_tent = (1./dt)*inner(u_tent - u0, v)*dx + inner(dot(u0.dx(0), u0), v)*dx + nu*inner((u_tent.dx(0) + u0.dx(0)), v.dx(0)) * dx + inner(p0.dx(0), v) * dx\
    - inner(f, v)*dx     # solve to u_
a_tent, L_tent = system(F_tent)
# step 2
u_tent_computed = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
F_p = inner(grad(p-p0), grad(q))*dx + (1./dt)*u_tent_computed.dx(0)*q*dx # + 2*p.dx(0)*q*ds(1) # tried to force dp/dn=0 on inflow
# TEST: prescribe Neumann outflow BC
# F_p = inner(grad(p-p0), grad(q))*dx + (1./dt)*u_tent_computed.dx(0)*q*dx + (1./dt)*(v_in_expr-u_tent_computed)*q*ds(2)
a_p, L_p = system(F_p)
A_p = assemble(a_p)
as_backend_type(A_p).set_nullspace(null_space)
print(A_p.array())
# step 2 rotation
# F_p_rot = inner(grad(p), grad(q))*dx + (1./dt)*u_tent_computed.dx(0)*q*dx + (1./dt)*(v_in_expr-u_tent_computed)*q*ds(2)

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import dx [as 别名]
    def step(self,
             dt,
             u1, p1,
             u, p0,
             u_bcs, p_bcs,
             f0=None, f1=None,
             verbose=True,
             tol=1.0e-10
             ):
        '''General pressure projection scheme as described in section 3.4 of
        :cite:`GMS06`.
        '''
        # Some initial sanity checkups.
        assert dt > 0.0

        # Define coefficients
        k = Constant(dt)

        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        # Compute tentative velocity step.
        # rho (u[-1] + (u.\nabla)u) = mu 1/r \div(r \nabla u) + rho g.
        with Message('Computing tentative velocity'):
            solver = NewtonSolver()
            solver.parameters['maximum_iterations'] = 5
            solver.parameters['absolute_tolerance'] = tol
            solver.parameters['relative_tolerance'] = 0.0
            solver.parameters['report'] = True
            # The nonlinear term makes the problem generally nonsymmetric.
            solver.parameters['linear_solver'] = 'gmres'
            # If the nonsymmetry is too strong, e.g., if u[-1] is large, then
            # AMG preconditioning might not work very well.
            # Use HYPRE-Euclid instead of ILU for parallel computation.
            solver.parameters['preconditioner'] = 'hypre_euclid'
            solver.parameters['krylov_solver']['relative_tolerance'] = tol
            solver.parameters['krylov_solver']['absolute_tolerance'] = 0.0
            solver.parameters['krylov_solver']['maximum_iterations'] = 1000
            solver.parameters['krylov_solver']['monitor_convergence'] = verbose

            ui = Function(self.W)
            step_problem = \
                TentativeVelocityProblem(ui,
                                         self.theta,
                                         self.rho, self.mu,
                                         u, p0, k,
                                         u_bcs,
                                         f0, f1,
                                         stabilization=self.stabilization,
                                         dx=self.dx
                                         )

            # Take u[-1] as initial guess.
            ui.assign(u[-1])
            solver.solve(step_problem, ui.vector())
            #div_u = 1/r * div(r*ui)
            #plot(div_u)
            #interactive()
        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        with Message('Computing pressure correction'):
            # The pressure correction is based on the update formula
            #
            #                           (    dphi/dr    )
            #     rho/dt (u_{n+1}-u*) + (    dphi/dz    ) = 0
            #                           (1/r dphi/dtheta)
            #
            # with
            #
            #     phi = p_{n+1} - p*
            #
            # and
            #
            #      1/r d/dr    (r ur_{n+1})
            #    +     d/dz    (  uz_{n+1})
            #    + 1/r d/dtheta(  utheta_{n+1}) = 0
            #
            # With the assumption that u does not change in the direction
            # theta, one derives
            #
            #  - 1/r * div(r * \nabla phi) = 1/r * rho/dt div(r*(u_{n+1} - u*))
            #  - 1/r * n. (r * \nabla phi) = 1/r * rho/dt  n.(r*(u_{n+1} - u*))
            #
            # In its weak form, this is
            #
            #   \int r * \grad(phi).\grad(q) *2*pi =
            #       - rho/dt \int div(r*u*) q *2*p
            #       - rho/dt \int_Gamma n.(r*(u_{n+1}-u*)) q *2*pi.
            #
            # (The terms 1/r cancel with the volume elements 2*pi*r.)
            # If Dirichlet boundary conditions are applied to both u* and u_n
            # (the latter in the final step), the boundary integral vanishes.
            #
            alpha = 1.0
            #alpha = 3.0 / 2.0
            rotational_form = False
            self._pressure_poisson(p1, p0,
                                   self.mu, ui,
                                   self.rho * alpha * ui / k,
                                   p_bcs=p_bcs,
                                   rotational_form=rotational_form,
                                   tol=tol,
                                   verbose=verbose
#.........这里部分代码省略.........