How to update the variables in the variational form object

[bw...@mst.edu]

I have a time-dependent problem to solve. So I have to update the solution after each time-step. If I define the variational form outside the loop, it seems that it won't get updated.  But if I write the form inside the loop, it seems to be very time-consuming and the computational cost is huge. Finally it'll stop somewhere before the time loop ends. I think this is due to run-out of memory. Does anyone have any good suggestions? Thanks in advance.  Here are the codes:

#phase field approach for 2D fracture problem

# A bar fixed on the left end and the right end is subject to boundary condition u=t*l

from dolfin import *

from mshr import *

from fenics import *

import numpy

# Geometry

l= 1.0; h = 0.1

# Material constant

E, nu = Constant(10.0E3), Constant(0.3)

#Lame constants

lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))

mu = Constant(E/(2*(1+nu)))

bulk = Constant(3*lmbda+2*mu);

k=Constant(bulk/3) #bulk modulus

gc = 0.25*1.0;

Gc = Constant(gc)

lc=Constant(0.1)

k_ell = Constant(1.e-6) # residual stiffness

# Loading

ut = 1. # reference value for the loading (imposed displacement)

f = 0. # bulk load

t=0.0

H_n=Constant(0.0)#previous history variable

 # Strain and stress

def eps(u):

    return sym(grad(u))

def tr_e_plus(u,un): #extensive/positive principal components of the strain tensor

    t=tr(eps(un))

    s= conditional(gt(t,0.0),tr(eps(u)),0.0)

    return s

def tr_e_minus(u,un): #extensive/positive principal components of the strain tensor

    t=tr(eps(un))

    s= conditional(lt(t,0.0),tr(eps(u)),0.0)

    return s

def g(d):  #stress degradation function

    return (1-d)**2*(1-k_ell)+k_ell

def sigma_t(u,dn,un):

    str_ele = 0.5*(nabla_grad(u) + nabla_grad(u).T)

    ICp = tr_e_plus(u,un)

    ICm = tr_e_minus(u,un)

    return k*(ICm  + g(dn)*ICp  )*Identity(2)+ g(dn)*2*mu*(str_ele -1/3*(ICm+ICp)*Identity(2))

def en_dens_p(u):#positive part of energy density

    str_ele =eps(u)

    IC = tr(str_ele)

    ICC = tr(str_ele * str_ele)

    I=conditional(gt(IC,0.0),IC,0.0)

    return (bulk/6*I**2) + mu*(ICC-1/3*IC**2)

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

# Define boundary sets for boundary conditions

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

def left(x, on_boundary):

    return near(x[0] * 0.01, 0)

        

def right(x, on_boundary):

    return near((x[0] - 1) * 0.01, 0.)

       

#==================================

# Create mesh 

#===================================

nx =100

ny =10

#Create a uniform finite element mesh of size 100*10

mesh = RectangleMesh(Point(0, 0), Point(l, h), nx, ny)

#==================================

# Define function space

#===================================

V = FunctionSpace(mesh, 'CG', 1)

W=VectorFunctionSpace(mesh, "CG", 1)

u , v = TrialFunction(W), TestFunction(W)

d, e=TrialFunction(V), TestFunction(V)

# Dirichlet boundary condition for a traction test boundary

zero_v=Constant((0.0,0.0))

u_L = zero_v

u_R = Expression(("t","0."),t = 0,degree=1)

# bc - u (imposed displacement)

Gamma_u_0 = DirichletBC(W, u_L, left)

Gamma_u_1 = DirichletBC(W, u_R, right)

bc_u = [Gamma_u_0, Gamma_u_1]

# bc - d (zero damage)

Gamma_d_0 = DirichletBC(V, 0.0, left)

Gamma_d_1 = DirichletBC(V, 0.0, right)

bc_d = [Gamma_d_0, Gamma_d_1]

d0=Constant(0.0)  #damage 

d_n = interpolate(d0, V)

u0=zero_v

u_n = interpolate(u0, W)

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

# Define variational problem

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

# Create VTK file for saving solution

vtkfile1=File('results/displacement.pvd')

vtkfile2=File('results/damage.pvd')

a_u=inner(sigma_t(u,d_n,u_n),grad(v))*dx

L_u=inner(Constant(0.0)*Identity(2),grad(v))*dx

a_d=(2*H_n+gc/lc)*d*e*dx+gc*lc*dot(grad(d), grad(e))*dx

L_d=2*H_n*e*dx

p = Function(W)#displacements

a=Function(V)#damage

num_steps=50

dt=1.0/num_steps

for i in range(num_steps):# Update current time

    t = t+dt

    Au= assemble(a_u)

    bu= assemble(L_u)

    solve(Au, p.vector(), bu)

    u_n.assign(p)

    H=conditional(gt(en_dens_p(u_n),H_n),en_dens_p(u_n),H_n)

    H_n=H

    #history variable

    print t

    Ad= assemble(a_d)

    bd= assemble(L_d)

    solve(Ad, a.vector(), bd)

    d_n.assign(a)

    u_R.t=t

    plot(p)

    plot(a)

    vtkfile1<< (p, t)

    vtkfile2<< (a, t)

[Carl Lundholm]

Yes, have a look at the demo for the advection-diffusion equation.

/Carl

---

Från: fenics-...@googlegroups.com <fenics-...@googlegroups.com> för bw...@mst.edu <bw...@mst.edu>
Skickat: den 22 maj 2018 20:36:10
Till: fenics-support
Ämne: [fenics-support] How to update the variables in the variational form object

 

--
You received this message because you are subscribed to the Google Groups "fenics-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fenics-suppor...@googlegroups.com.
To post to this group, send email to fenics-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fenics-support/c014d16c-2214-401a-a1d4-c32ae1b1e665%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Carl Lundholm]

---

Från: Carl Lundholm
Skickat: den 23 maj 2018 10:23
Till: Bo Wang
Ämne: SV: [fenics-support] How to update the variables in the variational form object

 

If you have to update something in every time-step then there's no way around it that I know of. 

But perhaps you could see if there are some things that you can move outside the time-stepping loop. There are usually more alternatives than "have all forms outside the time-stepping loop" OR "have all forms inside the time-stepping loop".

Some tips:

Having an assembly of a matrix at every time-step can be quite costly. Ideally you would like to just assemble the matrix once before the time-stepping and reuse the factorization as you've seen in the advection-diffusion demo. If the time-dependency in the matrix just concerns a spatially invariant parameter then the update should just concern a scalar-matrix multiplication. You should try to move out as much as possible from the time-stepping, e.g., your L_u is just constantly = 0 since it contains a factor Constant(0.0) so there's no need to assemble it in the time-stepping. Perhaps there's more? 

/Carl

---

Från: Bo Wang <bw...@mst.edu>
Skickat: den 22 maj 2018 22:32:06
Till: Carl Lundholm
Ämne: Re: [fenics-support] How to update the variables in the variational form object

 

Thanks. But my problem is if there is a simpler way to update the variables in the bi-linear  form (a) and linear form  (l). 

On Tue, May 22, 2018 at 2:44 PM, Carl Lundholm <car...@chalmers.se> wrote:

Yes, have a look at the demo for the advection-diffusion equation.

/Carl

---

Från: fenics-support@googlegroups.com <fenics-support@googlegroups.com> för bw...@mst.edu <bw...@mst.edu>

To unsubscribe from this group and stop receiving emails from it, send an email to fenics-support+unsubscribe@googlegroups.com.
To post to this group, send email to fenics-support@googlegroups.com.

[Carl Lundholm]

No problem and take a look at the the advection diffusion demo if you haven't. There they have the costly assembly of the system matrix outside the time-stepping. Maybe that can help.

/Carl

---

Från: Bo Wang <bw...@mst.edu>
Skickat: den 24 maj 2018 05:01:07
Till: Carl Lundholm
Ämne: Re: [fenics-support] How to update the variables in the variational form object

 

Thanks for your reply. Actually I'm quite new to fenics. But I think you are right. I'll try to use other methods like LU factorization to give it a try.And your tips helped a lot. Thanks very much.

On Wed, May 23, 2018 at 4:23 AM, Carl Lundholm <car...@chalmers.se> wrote:

If you have to update something in every time-step then there's no way around it that I know of. 

But perhaps you could see if there are some things that you can move outside the time-stepping loop. There are usually more alternatives than "have all forms outside the time-stepping loop" OR "have all forms inside the time-stepping loop".

Some tips:

Having an assembly of a matrix at every time-step can be quite costly. Ideally you would like to just assemble the matrix once before the time-stepping and reuse the factorization as you've seen in the advection-diffusion demo. If the time-dependency in the matrix just concerns a spatially invariant parameter then the update should just concern a scalar-matrix multiplication. You should try to move out as much as possible from the time-stepping, e.g., your L_u is just constantly = 0 since it contains a factor Constant(0.0) so there's no need to assemble it in the time-stepping. Perhaps there's more? 

/Carl

---

Från: Bo Wang <bw...@mst.edu>
Skickat: den 22 maj 2018 22:32:06
Till: Carl Lundholm
Ämne: Re: [fenics-support] How to update the variables in the variational form object

 

Thanks. But my problem is if there is a simpler way to update the variables in the bi-linear  form (a) and linear form  (l). 

On Tue, May 22, 2018 at 2:44 PM, Carl Lundholm <car...@chalmers.se> wrote:

Yes, have a look at the demo for the advection-diffusion equation.

/Carl

---

Från: fenics-...@googlegroups.com <fenics-support@googlegroups.com> för bw...@mst.edu <bw...@mst.edu>

To unsubscribe from this group and stop receiving emails from it, send an email to fenics-support+unsubscribe@googlegroups.com.