https://github.com/hpc-maths/GenEO
Revision 3881e864d0115355ab0398e22b489d3d1df7997a authored by gouarin on 30 March 2018, 12:29:37 UTC, committed by gouarin on 30 March 2018, 12:29:37 UTC
1 parent eb1bfe3
Tip revision: 3881e864d0115355ab0398e22b489d3d1df7997a authored by gouarin on 30 March 2018, 12:29:37 UTC
fix mpcg
fix mpcg
Tip revision: 3881e86
test_cg.py
from __future__ import print_function, division
import sys, petsc4py
petsc4py.init(sys.argv)
import mpi4py.MPI as mpi
from petsc4py import PETSc
import numpy as np
from elasticity import *
def rhs(coords, rhs):
rhs[..., 1] = -9.81
OptDB = PETSc.Options()
Lx = OptDB.getInt('Lx', 10)
Ly = OptDB.getInt('Ly', 1)
n = OptDB.getInt('n', 16)
nx = OptDB.getInt('nx', Lx*n)
ny = OptDB.getInt('ny', Ly*n)
hx = Lx/(nx - 1)
hy = Ly/(ny - 1)
da = PETSc.DMDA().create([nx, ny], dof=2, stencil_width=1)
da.setUniformCoordinates(xmax=Lx, ymax=Ly)
# constant young modulus
E = 30000
# constant Poisson coefficient
nu = 0.4
lamb = (nu*E)/((1+nu)*(1-2*nu))
mu = .5*E/(1+nu)
x = da.createGlobalVec()
b = buildRHS(da, [hx, hy], rhs)
A = buildElasticityMatrix(da, [hx, hy], lamb, mu)
A.assemble()
bcApplyWest(da, A, b)
ksp = PETSc.KSP().create()
ksp.setOperators(A)
ksp.setType("cg")
ksp.setInitialGuessNonzero(True)
pc = ksp.pc
pc.setType('none')
ksp.setFromOptions()
ksp.solve(b, x)
cg(A, b)
Computing file changes ...
