search

underworldcode / underworld2

underworld2: A parallel, particle-in-cell, finite element code for Geodynamics.
http://www.underworldcode.org/
Other
168 stars 58 forks source link

Neumann BC in visco-elastic problem #558

Open Yidali26 opened 3 years ago

Yidali26 commented 3 years ago

Hi there, I tried to test the possibility of Neumann Boundary condition when the material is visco-elastic. The test I set up is a simple shear problem with eta=1, miu=1in a 1x1 box. Left and right BC are periodic, bottom is no slip and top boundary is tau_xz=t. This configuration should give a velocity field on the top to be vx=1+t, but the model return some weird results. Below is the code:

import underworld as uw
from underworld import function as fn
import numpy as np
import matplotlib.pyplot as plt
from time import time as timer
elementType = "Q1/dQ0"
resX =32 #400
resY = 2
resZ = 32 #230
maxX = 1.
maxY = 0.02
maxZ = 1.

mesh = uw.mesh.FeMesh_Cartesian( elementType = (elementType),
                                 elementRes  = (resX,resY, resZ),
                                 minCoord    = (0., 0., 0.),
                                 maxCoord    = (maxX, maxY, maxZ),
                                 periodic    = [True, False, False]  )

velocityField    = mesh.add_variable(         nodeDofCount=mesh.dim )
appliedTractionField = uw.mesh.MeshVariable( mesh=mesh,    nodeDofCount=mesh.dim )
pressureField    = mesh.subMesh.add_variable( nodeDofCount=1 )
velocityField.data[:] = [0., 0., 0.]
pressureField.data[:] = 0.
swarm         = uw.swarm.Swarm(mesh, particleEscape=True)
swarmLayout   = uw.swarm.layouts.PerCellSpaceFillerLayout( swarm=swarm,particlesPerCell=16 )
pop_control = uw.swarm.PopulationControl(swarm, aggressive=True, particlesPerCell=16)
swarm.populate_using_layout( layout=swarmLayout )
previousStress         = swarm.add_variable( dataType="double", count=6 )
previousStress.data[:] = [0., 0., 0.,0., 0., 0.]
dStress         = swarm.add_variable( dataType="double", count=6 )
dStress.data[:] = [0., 0., 0.,0., 0., 0.]
advector        = uw.systems.SwarmAdvector( swarm=swarm,            velocityField=velocityField, order=2 )

iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
kWalls = mesh.specialSets["MinK_VertexSet"] + mesh.specialSets["MaxK_VertexSet"]
base = mesh.specialSets["MinK_VertexSet"]
top = mesh.specialSets["MaxK_VertexSet"]

vBC = uw.conditions.DirichletCondition( variable = velocityField, 
                                        indexSetsPerDof = (base, jWalls+kWalls, kWalls) )

appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data] = (0.,0.,0.)
nBC = uw.conditions.NeumannCondition( fn_flux=appliedTractionField, 
                                      variable=velocityField,
                                      indexSetsPerDof=(top, None, None) )

eta=1.
miu=1.
alpha = eta/miu
dt_e = eta/miu/10.

eta_eff = ( eta*dt_e ) / (alpha + dt_e)

strainRate = fn.tensor.symmetric( velocityField.fn_gradient )
tauHistoryFn = previousStress*eta_eff/(miu*dt_e)
viscoelasticeStressFn  =  eta_eff* (fn.misc.constant(2.) *strainRate+ previousStress/(miu*dt_e))

stokes = uw.systems.Stokes( velocityField = velocityField,
                               pressureField = pressureField,
                               voronoi_swarm = swarm,
                               conditions    = [vBC,nBC],
                               fn_viscosity  = eta_eff,
                               fn_stresshistory = tauHistoryFn
                          )
solver = uw.systems.Solver( stokes )
solver.set_inner_method(solve_type='lu')

def update():
    # Retrieve the maximum possible timestep for the advection system.
    dt = advector.get_max_dt()
    if dt > ( dt_e / 3. ):
        dt = dt_e / 3.
    advector.integrate(dt)
    newtime = time + dt

    # particle population control
    pop_control.repopulate()

    phi = dt / dt_e
    stressMapFn_data = viscoelasticeStressFn.evaluate(swarm)
    dStress.data[:] = (stressMapFn_data[:] - previousStress.data[:])/dt_e/2./miu
    previousStress.data[:] = ( phi*stressMapFn_data[:] + ( 1.-phi )*previousStress.data[:] )
    previousStress.data[:] = viscoelasticeStressFn.evaluate(swarm)

    return newtime, step+1

time = 0.
step = 0
t1=timer()
nsteps = 10
trec=np.zeros(nsteps)
vrec=np.zeros(nsteps)
sigmarec=np.zeros(nsteps)
while step<nsteps:
    appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time-eta_eff*previousStress.data[mesh.specialSets["MaxK_VertexSet"].data,4]/(miu*dt_e)
#     nBC = uw.conditions.NeumannCondition( fn_flux=appliedTractionField, 
#                                           variable=velocityField,
#                                           indexSetsPerDof=(top, None, None) )
#     stokes = uw.systems.Stokes( velocityField = velocityField,
#                                pressureField = pressureField,
#                                voronoi_swarm = swarm,
#                                conditions    = [vBC,nBC],
#                                fn_viscosity  = eta#_eff,
#                                #fn_stresshistory = tauHistoryFn
#                           )
#     solver = uw.systems.Solver( stokes )
    solver.solve()
    trec[step] = time
    vrec[step] = velocityField.data[velocityField.data.shape[0]-1,0]
    sigmarec[step] = previousStress.data[0,4]
    time, step = update()
    print(step,(timer()-t1)/60, time)

When configuring the boundary condition, I don't know if UW has already considered the elasticity term in Neumann condition, so I tried both , but neither works. The result of the code above is image

Yidali26 commented 3 years ago

If the elasticity is removed, the result agree with the analytical solution. I wonder if the elasticity could be compatible with the Neumann BC. Thanks

Yidali26 commented 3 years ago

Hi all, Is there any update or suggestion on this issue? I hope to set up some kind of open boundary with a visco-elastic problem, and this test is quite important to verify it's possible to set up the open boundary with a Neumann boundary. Thanks

lmoresi commented 3 years ago

There is something of a difficulty with non-linear stress boundary conditions / stress boundary conditions in which the element rheology plays a role because the usual integration by parts to create a surface integral in the weak form is not strictly appropriate.

However, what you are showing looks more like the elastic stresses are not being considered in the overall stress balance in the boundary condition - an inconsistency.

Any thoughts from @julesghub ?

Prof Louis Moresi

@.***

www.moresi.infohttp://www.moresi.info/

www.geo-down-under.org.auhttps://www.geo-down-under.org.au

@LouisMoresihttps://twitter.com/LouisMoresi

On 10 Aug 2021, 8:43 AM +1000, Yidali26 @.***>, wrote:

Hi all, Is there any update or suggestion on this issue? I hope to set up some kind of open boundary with a visco-elastic problem, and this test is quite important to verify it's possible to set up the open boundary with a Neumann boundary. Thanks

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/underworld2/issues/558#issuecomment-895599212, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADABPI2GK36IYPCGICE73Q3T4BKYNANCNFSM5BURHBVA.

Yidali26 commented 3 years ago

There is something of a difficulty with non-linear stress boundary conditions / stress boundary conditions in which the element rheology plays a role because the usual integration by parts to create a surface integral in the weak form is not strictly appropriate. However, what you are showing looks more like the elastic stresses are not being considered in the overall stress balance in the boundary condition - an inconsistency. Any thoughts from @julesghub ? Prof Louis Moresi @. www.moresi.info<http://www.moresi.info/> www.geo-down-under.org.au<https://www.geo-down-under.org.au> @LouisMoresihttps://twitter.com/LouisMoresi On 10 Aug 2021, 8:43 AM +1000, Yidali26 @.>, wrote: Hi all, Is there any update or suggestion on this issue? I hope to set up some kind of open boundary with a visco-elastic problem, and this test is quite important to verify it's possible to set up the open boundary with a Neumann boundary. Thanks — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#558 (comment)>, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADABPI2GK36IYPCGICE73Q3T4BKYNANCNFSM5BURHBVA.

Hi Louis, Thanks for your reply. Would the elasticity introduce any kind of non-linearity to this problem? I checked the weak form shown below, seems like no special treatment is required for the boundary condition with introduced elasticity, as the constitutive relation is not used in the boundary integration: webwxgetmsgimg This means my previous calculation that consider the tauhistory on the boundary traction must be wrong. Below is the new result with <img src="https://render.githubusercontent.com/render/math?math=\sigma{xz}=t"> and , which should return the velocitity and . The stress is correct while the velocity is a bit off from the analytical model.

Yidali26 commented 3 years ago

t_1sigma t_1vx t_2sigma t_2vx

Yidali26 commented 3 years ago

In the code appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time-eta_effpreviousStress.data[mesh.specialSets["MaxK_VertexSet"].data,4]/(miudt_e) should be replaced with appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time

julesghub commented 3 years ago

Hi @Yidali26, sorry for the late reply.

In the code appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time-eta_eff_previousStress.data[mesh.specialSets["MaxK_VertexSet"].data,4]/(miu_dt_e) should be replaced with appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time

You're right! There is a data structure inconsistency in this line. previousStress.data[mesh.specialSets["MaxK_VertexSet"].data,0]

previousStress.data[] is a swarm variables and is referenced by particle indices. mesh.specialSets["MaxK_VertextSet"] are the indices of the nodes, not particles. You can't directly use mesh.specialSets[] to index swarm variables. This was most likely the issue.

Also I see the model is 3D but only 2 elements wide in one direction. Perhaps remove the thin direction for testing.

Yidali26 commented 3 years ago

Hi @Yidali26, sorry for the late reply.

In the code appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time-eta_eff_previousStress.data[mesh.specialSets["MaxK_VertexSet"].data,4]/(miu_dt_e) should be replaced with appliedTractionField.data[mesh.specialSets["MaxK_VertexSet"].data,0] = time

You're right! There is a data structure inconsistency in this line. previousStress.data[mesh.specialSets["MaxK_VertexSet"].data,0]

previousStress.data[] is a swarm variables and is referenced by particle indices. mesh.specialSets["MaxK_VertextSet"] are the indices of the nodes, not particles. You can't directly use mesh.specialSets[] to index swarm variables. This was most likely the issue.

Also I see the model is 3D but only 2 elements wide in one direction. Perhaps remove the thin direction for testing.

Hi @julesghub, Thanks for your reply. Yes I set up a very thin layer in y dimension because I do the test for a 3D problem. Below is the result from a similar 2D problem, where I just removed the y dimension in the previous code. Also need to mention in my previous post I had a redundant 1/2 factor in front of the analytical vx. So the vx=t+1 for sigma=t and vx=t^2+2t for sigma =t^2. The results still looks like the previous ones, that the sigma is correct but vx is incorrect. t_1sigma t_1vx

t_2sigma t_2vx

code: 

resX =32 #400
resY = 2
resZ = 32 #230
maxX = 1.
maxY = 0.02
maxZ = 1.

mesh = uw.mesh.FeMesh_Cartesian( elementType = (elementType),
                                 elementRes  = (resX, resZ),
                                 minCoord    = (0.,  0.),
                                 maxCoord    = (maxX,  maxZ),
                                 periodic    = [True,  False]  )

velocityField    = mesh.add_variable(         nodeDofCount=mesh.dim )
appliedTractionField = uw.mesh.MeshVariable( mesh=mesh,    nodeDofCount=mesh.dim )
pressureField    = mesh.subMesh.add_variable( nodeDofCount=1 )
velocityField.data[:] = [0., 0.]
pressureField.data[:] = 0.
swarm         = uw.swarm.Swarm(mesh, particleEscape=True)
swarmLayout   = uw.swarm.layouts.PerCellSpaceFillerLayout( swarm=swarm,particlesPerCell=16 )
pop_control = uw.swarm.PopulationControl(swarm, aggressive=True, particlesPerCell=16)
swarm.populate_using_layout( layout=swarmLayout )
previousStress         = swarm.add_variable( dataType="double", count=3 )
previousStress.data[:] = [0., 0., 0.]
dStress         = swarm.add_variable( dataType="double", count=3 )
dStress.data[:] = [0., 0., 0.]
advector        = uw.systems.SwarmAdvector( swarm=swarm,            velocityField=velocityField, order=2 )

iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
base = mesh.specialSets["MinJ_VertexSet"]
top = mesh.specialSets["MaxJ_VertexSet"]

vBC = uw.conditions.DirichletCondition( variable = velocityField, 
                                        indexSetsPerDof = (base, jWalls) )

appliedTractionField.data[mesh.specialSets["MaxJ_VertexSet"].data] = (0.,0.)
nBC = uw.conditions.NeumannCondition( fn_flux=appliedTractionField, 
                                      variable=velocityField,
                                      indexSetsPerDof=(top, None) )

eta=1.
miu=1.
alpha = eta/miu
dt_e = eta/miu/10.

eta_eff = ( eta*dt_e ) / (alpha + dt_e)

strainRate = fn.tensor.symmetric( velocityField.fn_gradient )
tauHistoryFn = previousStress*eta_eff/(miu*dt_e)
viscoelasticeStressFn  =  eta_eff* (fn.misc.constant(2.) *strainRate+ previousStress/(miu*dt_e))

stokes = uw.systems.Stokes( velocityField = velocityField,
                               pressureField = pressureField,
                               voronoi_swarm = swarm,
                               conditions    = [vBC,nBC],
                               fn_viscosity  = eta_eff,
                               fn_stresshistory = tauHistoryFn
                          )
solver = uw.systems.Solver( stokes )
solver.set_inner_method(solve_type='lu')

def update():
    # Retrieve the maximum possible timestep for the advection system.
    dt = advector.get_max_dt()
    if dt > ( dt_e / 3. ):
        dt = dt_e / 3.
    advector.integrate(dt)
    newtime = time + dt

    # particle population control
    pop_control.repopulate()

    phi = dt / dt_e
    stressMapFn_data = viscoelasticeStressFn.evaluate(swarm)
    dStress.data[:] = (stressMapFn_data[:] - previousStress.data[:])/dt_e/2./miu
    previousStress.data[:] = ( phi*stressMapFn_data[:] + ( 1.-phi )*previousStress.data[:] )
    previousStress.data[:] = viscoelasticeStressFn.evaluate(swarm)

    return newtime, step+1

time = 0.
step = 0
t1=timer()
nsteps = 50
trec=np.zeros(nsteps)
vrec=np.zeros(nsteps)
sigmarec=np.zeros(nsteps)
while step<nsteps:
    appliedTractionField.data[mesh.specialSets["MaxJ_VertexSet"].data,0] = time**2#-eta_eff*previousStress.data[mesh.specialSets["MaxK_VertexSet"].data,4]/(miu*dt_e)
#     nBC = uw.conditions.NeumannCondition( fn_flux=appliedTractionField, 
#                                           variable=velocityField,
#                                           indexSetsPerDof=(top, None, None) )
#     stokes = uw.systems.Stokes( velocityField = velocityField,
#                                pressureField = pressureField,
#                                voronoi_swarm = swarm,
#                                conditions    = [vBC,nBC],
#                                fn_viscosity  = eta_eff,
#                                fn_stresshistory = tauHistoryFn
#                           )
#     solver = uw.systems.Solver( stokes )
    solver.solve()
    trec[step] = time
    vrec[step] = velocityField.data[velocityField.data.shape[0]-1,0]
    sigmarec[step] = previousStress.data[0,2]
    time, step = update()
    print(step,(timer()-t1)/60, time)

Do you have any idea about why the vx doesn't give the correct solution? Looking forward to your reply. Thanks Yida

julesghub commented 3 years ago

Hi @Yidali26, Thanks for the 2D version. Some thoughts on this:

  1. The initial 1st solve is solving for 0 external force, so our stokes solvers will struggle to find a suitable solution. The viscoelastic stresses (solution dependent) may be spurious and could be causing the oscillations for the vx solution.
  2. You can't define dirichlet conditions on nodes with periodic velocity conditions, i.e. on the iWalls. So define base to NOT include the periodic iWall nodes. ie. base = mesh.specialSets["MinJ_VertexSet"] - iWall. maybe try the same for the top too? Something to test.
Yidali26 commented 3 years ago

Hi @Yidali26, Thanks for the 2D version. Some thoughts on this:

1. The initial 1st solve is solving for 0 external force, so our stokes solvers will struggle to find a suitable solution. The viscoelastic stresses (solution dependent) may be spurious and could be causing the oscillations for the vx solution.

2. You can't define dirichlet conditions on nodes with periodic velocity conditions, i.e. on the iWalls. So define `base` to NOT include the periodic iWall nodes. ie.
   `base = mesh.specialSets["MinJ_VertexSet"] - iWall`.
   maybe try the same for the top too? Something to test.

Hi Julian, Thanks for your suggestions. For point 1, I tested the code and the code just return the v=[0.,0.] without any problem. For this set up, at t=0 the top bc is \tau_xz=0, but not a ill-constraint problem. As both internal forces and boundary tractions are zeros, the RHS term in the weak form might be a zero vector, but I'm using a LU solver, so I suppose there wouldn't be any convergence issue. For point 2. I changed the code according to your suggestion, but seem like the problem still exist. Below is the results, quite similar to the earlier ones. image image Looking forward to your reply! Thanks Yida