Updating time dependent expression leading to wrong output files

dolfinx
1

Hello Everyone,

I want to have a linearly increasing load on one of the edges and a constant value after time t = 20s. As long as I have constant load throughout, everything is fine. Bur when I update the load on the boundary with time (lines 108 - 119), the output file behaves weird. The time in paraview starts at my end time = 2000s and goes to 20001. Although there should be 200 steps. I also printed out the values in the time loop. I see that the time and the load is getting updated. However, I do not see why the output files are behaving weird. I am using Ubuntu 22 and dolfinx 0.7.2 and Paraview 5.12(latest). Any help/suggestions to fix this issue? I am attaching the minimum example below

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

from dolfinx.io import VTXWriter
from dolfinx import nls
import numpy as np
import ufl
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

from mpi4py import MPI
from dolfinx import fem, mesh, plot
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc #solvers

#-------------------------------------------------------------------#

# Geometry
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [10.0, 10.0]], [20, 20], mesh.CellType.quadrilateral)

#-------------------------------------------------------------------#

# Function Spaces and Mixed Space
u_el = ufl.VectorElement("CG", domain.ufl_cell(), 2) 
p_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) 
element = ufl.MixedElement([u_el, p_el]) 
V = fem.FunctionSpace(domain, element) 

up = fem. Function(V) 
u_sp = fem.FunctionSpace(domain, u_el) 
p_sp = fem.FunctionSpace(domain, p_el) 

Vu, up_to_u = V.sub(0).collapse() 
u_n = fem. Function(Vu) 
Vp, V_to_Vp = V.sub(1).collapse()

u, p = ufl.split(up) 
(del_u, del_p) = ufl.TestFunctions(V)

#-------------------------------------------------------------------#

# Define temporal parameters
t = 0  # Start time
T = 2000.0  # Final time
num_steps = 200
dt = T / num_steps  # time step size

#-------------------------------------------------------------------#

# Boundary Detection
boundaries = [
        (1, lambda x: np.isclose(x[0], 0)),
        (2, lambda x: np.isclose(x[1], 0)),
        (3, lambda x: np.isclose(x[0], 10.0)),
        (4, lambda x: np.isclose(x[1], 10.0)),]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
dim = domain.topology.dim

for (marker, locator) in boundaries:
    facets = mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Boundary Conditions

def bc_vec_zero(x):
    return np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)

def bc_scalar_zero(x):
    return np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)

list_bc = []

# Displacement Boundary
u_bc = fem.Function(u_sp)
u_bc.interpolate(bc_vec_zero)

facets = facet_tag.indices[facet_tag.values == 2]
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(1), u_sp.sub(1)), fdim, facets) 
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))

facets = facet_tag.indices[facet_tag.values == 1]
left_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets) 
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))

facets = facet_tag.indices[facet_tag.values == 3]
right_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets) 
list_bc.append(fem.dirichletbc(u_bc, right_dofs, V.sub(0)))

# Pressure Boundary
p_bc = fem.Function(p_sp)
p_bc.interpolate(bc_scalar_zero)

facets = facet_tag.indices[facet_tag.values == 4]
top_dofs = fem.locate_dofs_topological((V.sub(1), p_sp), fdim, facets) 
list_bc.append(fem.dirichletbc(p_bc, top_dofs, V.sub(1)))

# Load Boundary

def load_func(t):
    if t == 0:
        eval = 0
    elif t <= 20.0:
        eval = -(10000/20) * t
    else:
        eval = -10000
    return eval

load = fem.Constant(domain, ScalarType(load_func(t)))
print(load.value)
load_bc = load * ufl.FacetNormal(domain)

#-------------------------------------------------------------------#

# Parameters

# Material Parameters
mat_k_FOS = 1e-4
mat_gamma_FR = 1e4
mat_rho_SR = 2e3
mat_rho_FR = 1e3
mat_n_SOS = 0.67

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
b = ufl.variable(F * F.T)

# Invariants of deformation tensors
J = ufl.variable(ufl.det(F))

# Solid Velocity
vel_s = (u - u_n) / dt

# Darcy
dar = - (ufl.grad(p))

# Stress
mu_1 = default_scalar_type(5.58e6)
lmbda_1 = default_scalar_type(8.38e6)
P_ef_2 = (1/J) * (mu_1 * (b - I) + lmbda_1 * ufl.ln(J) * I)
P = P_ef_2 - p*I

# Integration measures
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)
dx = ufl.Measure("dx", domain=domain)

# Weak Form
Mom_B = ufl.inner(P, ufl.grad(del_u)) * ufl.dx - ufl.inner(del_u, load_bc) * ds(4)
Vol_B = (ufl.div(vel_s) * del_p + ufl.inner(dar, ufl.grad(del_p))) * ufl.dx

weak_form = Mom_B + Vol_B

#-------------------------------------------------------------------#

# Solver
# Initialise non-linear problem
problem = NonlinearProblem(weak_form, up, list_bc)

# Initialise newton solver
solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.max_it = 100000

# Configure mumps 
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Initialize history values
u_n.x.array[:] = 0.0

#-------------------------------------------------------------------#

# --- Solve problem ---
duration_solve = 0.0
for n in range(num_steps):
    # Update time
    t += dt

    # Update time
    load.value = load_func(t)
    print(load.value)

    # Solve newton-steps
    duration_solve -= MPI.Wtime()
    niter, converged = solver.solve(up)
    duration_solve += MPI.Wtime()

    PETSc.Sys.Print(
        "Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
            t, duration_solve, niter
        )
    )

    u_n.x.array[:] = up.x.array[up_to_u]

    u = up.sub(0).collapse()
    p = up.sub(1).collapse()

    u.name = "disp"
    file1 = dolfinx.io.VTXWriter(domain.comm, "output1.bp", u, "bp4")
    file1.write(t)
    print(t)
    file1.close()

    p.name = "pressure"
    file2 = dolfinx.io.VTXWriter(domain.comm, "output2.bp", p, "bp4")
    file2.write(t)
    file2.close()

Thank you in advance :slight_smile:

1 Like
2

You should move the definition of file1 and file2 outside the temporal loop, as you are currently overwriting data in an unclear fashion.
i.e.

file1 =  dolfinx.io.VTXWriter(domain.comm, "output1.bp", u, "bp4")
for n in range(num_steps):
    # Add code here

    file1.write(t)
# After loop
file1.close()

Please also try to simplify your problem (currently there is a lot of code to parse).

1 Like
3

Thank you very much for the quick response. I made the changes but I have the following error now regarding the first line itself.

file1 = dolfinx.io.VTXWriter(domain.comm, "output1.bp", u, "bp4")
file2 = dolfinx.io.VTXWriter(domain.comm, "output2.bp", p, "bp4")

for n in range(num_steps):
    # Update time
    t += dt

    # Update load
    load. Value = load_func(t)

    # Solve newton-steps
    duration_solve -= MPI.Wtime()
    niter, converged = solver.solve(up)
    duration_solve += MPI.Wtime()

    PETSc.Sys.Print(
        "Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
            t, duration_solve, niter
        )
    )

    u_n.x.array[:] = up.x.array[up_to_u]

    u = up.sub(0).collapse()
    p = up.sub(1).collapse()

    u.name = "disp"
    p.name = "pressure"

    file1.write(t)
    file2.write(t)

file1.close()
file2.close()

Exception has occurred: AttributeError
'Indexed' object has no attribute 'function_space'
AttributeError: 'ListTensor' object has no attribute 'geometry'

During handling of the above exception, another exception occurred:

AttributeError: 'ListTensor' object has no attribute 'function_space'

During handling of the above exception, another exception occurred:

  File "/home/igupta/studies/Basic_bi_copy/basic_bi.py", line 206, in <module>
    file1 = dolfinx.io.VTXWriter(domain.comm, "output1.bp", u, "bp4")
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'Indexed' object has no attribute 'function_space'

Thank you for the suggestion. I will reduce the code further for future questions.

4

Here

u and p are still the TrialFunctions, and not yet the solution.

One option would be to define a u_solution and p_solution as dolfinx.fem.Function (for your own sanity you may not want to call them again u and p), and then update them inside the loop as done e.g. in dolfinx/python/demo/demo_navier-stokes.py at v0.7.2 · FEniCS/dolfinx · GitHub

EDIT: I hadn’t realize that your code was using a MixedElement, and not block matrices/vectors. The link above doesn’t cover that case, but should still give you an idea about what to do.

5

Thanks for the link. I did the following. I created fem. Function using the sub().collapse() to update it in the time loop. But still, the output file is giving wrong results like in the picture below for pressure. What could be wrong here? (Sorry if it is a stupid question. I am quite new to Fenics) .

Vu_sol, up_to_u_sol = V.sub(0).collapse() 
u_sol = fem. Function(Vu_sol) 

Vp_sol, up_to_p_sol = V.sub(1).collapse() 
p_sol = fem. Function(Vp_sol) 

u_sol.name = "disp"
p_sol.name = "pressure"

file1 = dolfinx.io.VTXWriter(domain.comm, "output1.bp", u_sol, "bp4")
file2 = dolfinx.io.VTXWriter(domain.comm, "output2.bp", p_sol, "bp4")

# --- Solve problem ---
duration_solve = 0.0

for n in range(num_steps):
    # Update time
    t += dt

    # Update load
    load.value = load_func(t)
    print(load.value)
    print(t)

    # Solve newton-steps
    duration_solve -= MPI.Wtime()
    niter, converged = solver.solve(up)
    duration_solve += MPI.Wtime()

    PETSc.Sys.Print(
        "Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
            t, duration_solve, niter
        )
    )

    u_n.x.array[:] = up.x.array[up_to_u]

    u_eval = up. Sub(0).collapse()
    p_eval = up. Sub(1).collapse()

    u_sol.interpolate(u_eval)
    p_sol.interpolate(p_eval)

    file1.write(t)
    file2.write(t)

file1.close()
file2.close()

Screenshot from 2023-12-14 11-18-47
Screenshot from 2023-12-14 11-18-471066Ă—711 44.1 KB

6

It’s unclear to me what “wrong” means. The only blatantly wrong case is if the exported solution was the same the at all times, but the solution of the solve was indeed changing (as shown eg by printing its norm inside the time loop). Anything else is not a visualization issue, rather a problem with how you set up your problem.

Your code is very long and you haven’t given any detail on which kind of physical problem you are trying to solve, so it’s unlikely that anyone else apart form you will have an understanding of what you are trying to do and would be able to help.

1 Like
7

Thank you very much for the response. It was indeed an error with the weak formulation. On paper it was fine but I made a mistake while coding. Thank you for the advise earlier :slight_smile: