Hard bcs constraint with sigmoid approach

FEniCSx.py

+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, log
+# from dolfinx.fem import Measure
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import grad, dx, inner, TestFunction, as_vector
+from basix.ufl import element
+from petsc4py import PETSc
+import matplotlib.pyplot as plt
+from HeatCoefficients import *
+
+def solve_heatequation_dimensionless(T_ICE, T_MELTING, dt, t_end, T_REF, relax_param=None, number_cells=128, rtol=1e-4, max_iter=500):
+    x0 = 0
+    x1 = 1
+
+    # Define boundaries
+    def Left(x):
+        return np.isclose(x[0], x0)
+
+    def Right(x):
+        return np.isclose(x[0], x1)
+
+    # Define mesh
+    msh = mesh.create_interval(MPI.COMM_WORLD, number_cells, points=[x0,x1], dtype=np.float64)
+    
+    # Define Finite Element
+    CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+
+    # Define Mixed Function Space
+    W = fem.functionspace(msh, CG1_elem)
+
+    # Define Trial- and Testfunctions
+    Theta = fem.Function(W)
+    Theta_n = fem.Function(W)
+    v = TestFunction(W)
+
+    # Define boundary conditions
+    Theta_LEFT = 0.
+    Theta_RIGHT = 1.
+    def u_left_(x):
+        return np.full_like(x[0], Theta_LEFT)
+    def u_right_(x):
+        return np.full_like(x[0], Theta_RIGHT)
+
+    u_left_bcs = fem.Function(W)
+    u_left_bcs.interpolate(u_left_)
+    u_right_bcs = fem.Function(W)
+    u_right_bcs.interpolate(u_right_)
+
+    facet_left_dofs = fem.locate_dofs_geometrical(W, Left)
+    facet_right_dofs = fem.locate_dofs_geometrical(W, Right)
+    bc_left_T = fem.dirichletbc(PETSc.ScalarType(Theta_LEFT), facet_left_dofs, W)
+    bc_right_T = fem.dirichletbc(PETSc.ScalarType(Theta_RIGHT), facet_right_dofs, W)
+
+
+    # Define dimless properties
+    delta_T = T_MELTING - T_ICE  # ΔT
+    beta = delta_T / T_ICE  # β
+    rho_ref = rho_physical(T_REF) # ρ_ref
+    c_p_ref = c_p_physical(T_REF) # c_p_ref
+    lambda_ref = lambda_physical(T_REF) # λ_ref
+    L = x1 - x0
+    t_ref = (rho_ref * c_p_ref * L**2 ) / lambda_ref # t_ref
+
+    bcs = [bc_left_T, bc_right_T]
+
+    # Variational Formulation    
+    F = (
+        rho_star(Theta, delta_T=delta_T, T_REF=T_REF, rho_ref=rho_ref) 
+        * c_p_star(Theta, delta_T=delta_T, T_REF=T_REF, c_p_ref=c_p_ref) 
+        * (Theta - Theta_n) / dt * v * dx 
+        + lambda_star(Theta, delta_T=delta_T, T_REF=T_REF, lambda_ref=lambda_ref) 
+        * inner(grad(Theta), grad(v)) * dx
+    )
+    '''
+    F = a - L
+      = T * v * dx + dt * inner(grad(T), grad(v)) *dx - (T_n) * v * dx
+      = (T - T_n) * v * dx + dt * inner(grad(T), grad(v)) *dx
+      = (T - T_n) / dt * v * dx + inner(grad(T), grad(v)) *dx
+    '''
+
+    # Initialize initial guess for T to be 0°C = 273.15K
+    Theta_init_temp = Theta_LEFT
+    Theta_init = fem.Function(W)
+    Theta_init.interpolate(lambda x: np.full_like(x[0], Theta_init_temp))
+
+    with Theta.vector.localForm() as u_loc:
+        u_loc.set(0)
+    Theta.interpolate(Theta_init)
+    Theta_n.interpolate(Theta_init)
+
+    # Define Nonlinear Problem
+    problem = NonlinearProblem(F, Theta, bcs=bcs)
+
+    # Define Newton Solver
+    solver = NewtonSolver(MPI.COMM_WORLD, problem)
+    solver.convergence_criterion = "incremental"
+    solver.rtol = rtol
+    solver.relaxation_parameter = relax_param
+    solver.max_it = max_iter
+    solver.report = True
+
+    log.set_log_level(log.LogLevel.INFO)
+
+    cur_time = 0
+    # t_history, T_history = [cur_time], [np.full_like(T.vector, T_init_temp)]
+    x_grad = grad(Theta_n)[0]
+    expr = fem.Expression(x_grad, W.element.interpolation_points())
+    w = fem.Function(W)
+    w.interpolate(expr)
+    t_history, Theta_history, dTheta_history  = [cur_time], [np.array(Theta_init.vector)], [w.x.array]
+    while cur_time < t_end:
+        print(f'current time: {cur_time}')
+
+        n, converged = solver.solve(Theta)
+        assert (converged)
+
+        cur_time += dt
+        Theta.x.scatter_forward
+        Theta_n.x.array[:] = Theta.x.array
+
+        x_ = np.array(msh.geometry.x[:,0])
+        Theta_ = np.array(Theta_n.vector)
+
+        t_history.append(cur_time)
+        Theta_history.append(Theta_)
+        x_grad = grad(Theta_n)[0]
+        expr = fem.Expression(x_grad, W.element.interpolation_points())
+        w = fem.Function(W)
+        w.interpolate(expr)
+        dTheta_history.append(w.x.array)
+
+    x = msh.geometry.x[:,0]
+
+    return x, Theta_history, t_history, dTheta_history
+    
+    
+if __name__ == '__main__':
+    # Physical parameters
+    T_ICE = 100 # [K]
+    T_MELTING = 273.15 # [K]
+
+    T_REF = 100
+    t_end = 1
+    dt = .1
+
+    # Solver and mesh parameters
+    number_cells = 64
+    relax_param = 1.
+    return_type = 'Vector'
+    rtol = 1e-4
+    max_iter = 2_500
+
+    
+    x, Theta_history, t_history, dTheta_history = solve_heatequation_dimensionless(T_ICE, T_MELTING, dt, t_end, T_REF, relax_param, number_cells, rtol, max_iter)
+        
+    plt.figure()
+    for index, time in enumerate(t_history):
+        plt.plot(x, Theta_history[index])#, label=f'time = {time/60/60:.2f} hours')
+    plt.legend()
+    plt.grid()
+    plt.xlabel('x')
+    plt.ylabel('T')
+    plt.show()
+
+    # Visualize as a 3D plot
+    X, T = np.meshgrid(x, t_history)
+    Theta_history = np.array(Theta_history)
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    ax.plot_surface(X, T, Theta_history)
+    ax.set_xlabel('x')
+    ax.set_ylabel('t')
+    ax.set_zlabel('Theta')
+    plt.show()
+
+    # Visualize as a 3D plot
+    dTheta_history = np.array(dTheta_history)
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    ax.plot_surface(X, T, dTheta_history)
+    ax.set_xlabel('x')
+    ax.set_ylabel('t')
+    ax.set_zlabel('Theta')
+    plt.show()