Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
u2.create(grid, "", PISM.WITH_GHOSTS)
u2.copy_from(ssarun.ssa.solution())
Td_fd = PISM.IceModelVec2V()
Td_fd.create(grid, "", PISM.WITH_GHOSTS)
Td_fd.copy_from(u2)
Td_fd.add(-1, u1)
Td_fd.scale(1. / eps)
d_Td = PISM.IceModelVec2V()
d_Td.create(grid, "", PISM.WITH_GHOSTS)
d_Td.copy_from(Td_fd)
d_Td.add(-1, Td)
n_Td_fd = Td_fd.norm(PETSc.NormType.NORM_2)
n_Td_l2 = Td.norm(PETSc.NormType.NORM_2)
n_Td_l1 = Td.norm(PETSc.NormType.NORM_1)
n_Td_linf = Td.norm(PETSc.NormType.NORM_INFINITY)
n_d_Td_l2 = d_Td.norm(PETSc.NormType.NORM_2)
n_d_Td_l1 = d_Td.norm(PETSc.NormType.NORM_1)
n_d_Td_linf = d_Td.norm(PETSc.NormType.NORM_INFINITY)
PISM.verbPrintf(1, grid.com, "(i,j)=(%d,%d)\n" % (i, j))
PISM.verbPrintf(1, grid.com, "apply_linearization(d): l2 norm %.10g; finite difference %.10g\n" %
(n_Td_l2, n_Td_fd))
r_d_l2 = 0
if n_Td_l2 != 0:
r_d_l2 = n_d_Td_l2 / n_Td_l2
r_d_l1 = 0
if n_Td_l1 != 0:
drhs_fd.scale(1. / eps)
drhs = PISM.IceModelVec2V()
drhs.create(grid, "", PISM.WITHOUT_GHOSTS)
ssarun.ssa.apply_jacobian_design(u1, d, drhs)
d_drhs = PISM.IceModelVec2V()
d_drhs.create(grid, "", PISM.WITHOUT_GHOSTS)
d_drhs.copy_from(drhs)
d_drhs.add(-1, drhs_fd)
n_drhs_fd = drhs_fd.norm(PETSc.NormType.NORM_2)
n_drhs_l2 = drhs.norm(PETSc.NormType.NORM_2)
n_drhs_l1 = drhs.norm(PETSc.NormType.NORM_1)
n_drhs_linf = drhs.norm(PETSc.NormType.NORM_INFINITY)
n_d_drhs_l2 = d_drhs.norm(PETSc.NormType.NORM_2)
n_d_drhs_l1 = d_drhs.norm(PETSc.NormType.NORM_1)
n_d_drhs_linf = d_drhs.norm(PETSc.NormType.NORM_INFINITY)
PISM.verbPrintf(1, grid.com, "\nTest Jacobian Design (Comparison with finite differences):\n")
PISM.verbPrintf(1, grid.com, "jacobian_design(d): l2 norm %.10g; finite difference %.10g\n" %
(n_drhs_l2, n_drhs_fd))
if n_drhs_linf == 0:
PISM.verbPrintf(1, grid.com, "difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" %
(n_d_drhs_l2, n_d_drhs_l1, n_d_drhs_linf))
else:
PISM.verbPrintf(1, grid.com, "relative difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" %
(n_d_drhs_l2 / n_drhs_l2, n_d_drhs_l1 / n_drhs_l1, n_d_drhs_linf / n_drhs_linf))
view(d, d_viewer)
def iteration(self,it,eta,objVal,penaltyVal,d,diff_d,grad_d,u,diff_u,grad_u,grad):
print "----------------------------------------------------------"
print "Iteration %d" % it
print "RMS misfit: %g" % math.sqrt(penaltyVal)
print "sqrt(design objective) %g; weighted %g" % (math.sqrt(objVal),math.sqrt(objVal/eta))
print "gradient: design %g state %g sum %g " % (grad_d.norm(PETSc.NormType.NORM_2)/eta,grad_u.norm(PETSc.NormType.NORM_2),grad.norm(PETSc.NormType.NORM_2))
print "tikhonov functional: %g" % (objVal/eta + penaltyVal)
def iteration(self,it,eta,objVal,penaltyVal,d,diff_d,grad_d,u,diff_u,grad_u,grad):
print "----------------------------------------------------------"
print "Iteration %d" % it
print "RMS misfit: %g" % math.sqrt(penaltyVal)
print "sqrt(design objective) %g; weighted %g" % (math.sqrt(objVal),math.sqrt(objVal/eta))
print "gradient: design %g state %g sum %g " % (grad_d.norm(PETSc.NormType.NORM_2)/eta,grad_u.norm(PETSc.NormType.NORM_2),grad.norm(PETSc.NormType.NORM_2))
print "tikhonov functional: %g" % (objVal/eta + penaltyVal)
# print "doing iteration %d: gradients %g %g %g" % (it,grad_d.norm(PETSc.NormType.NORM_2),eta*grad_u.norm(PETSc.NormType.NORM_2),grad.norm(PETSc.NormType.NORM_2))
def norm(self, name):
if name == 'linf':
return self._core.norm(petsc4py.PETSc.NormType.NORM_INFINITY)
if name == 'l2':
return self._core.norm(petsc4py.PETSc.NormType.NORM_2)
if name == 'l1':
return self._core.norm(petsc4py.PETSc.NormType.NORM_1)
raise ValueError()
#
# x.assemblyBegin() # Needed in order to work on vector.
# x.assemblyEnd()
x = PETSc.Vec().createSeq(n) # Faster way to create a sequential vector.
x.setValues(range(n), range(n)) # x = [0 1 ... 9]
x.shift(1) # x = x + 1 (add 1 to all elements in x)
print 'Performing various vector operations on x =', x.getArray()
print 'Sum of elements of x =', x.sum()
print 'Dot product with itself =', x.dot(x)
print '1-norm =', x.norm(PETSc.NormType.NORM_1)
print '2-norm =', x.norm()
print 'Infinity-norm =', x.norm(PETSc.NormType.NORM_INFINITY)
print 'Minimum element in x (index, value) =', x.min()
print 'Maximum element in x (index, value) =', x.max()
def norm(self,name):
if name == 'linf':
return self._core.norm(PETSc.NormType.NORM_INFINITY)
if name == 'l2':
return self._core.norm(PETSc.NormType.NORM_2)
if name == 'l1':
return self._core.norm(PETSc.NormType.NORM_1)
raise ValueError()
y = da_fine.createGlobalVec()
# x_coarse.pointwiseMult(x_coarse)
# PETSc.Mat.Restrict(B, x_coarse, y)
B.mult(x_coarse, y)
# y.pointwiseMult(vec, y)
# PETSc.VecPointwiseMult()
# print(y.getArray())
# print(x_coarse.getArray())
print((y-x_fine).norm(PETSc.NormType.NORM_INFINITY))
y_coarse = da_coarse.createGlobalVec()
B.multTranspose(x_fine, y_coarse)
y_coarse.pointwiseMult(vec, y_coarse)
print((y_coarse - x_coarse).norm(PETSc.NormType.NORM_INFINITY))
# x.setSizes(n)
# x.setType('seq') # 'seq' means sequential vector.
#
# x.assemblyBegin() # Needed in order to work on vector.
# x.assemblyEnd()
x = PETSc.Vec().createSeq(n) # Faster way to create a sequential vector.
x.setValues(range(n), range(n)) # x = [0 1 ... 9]
x.shift(1) # x = x + 1 (add 1 to all elements in x)
print 'Performing various vector operations on x =', x.getArray()
print 'Sum of elements of x =', x.sum()
print 'Dot product with itself =', x.dot(x)
print '1-norm =', x.norm(PETSc.NormType.NORM_1)
print '2-norm =', x.norm()
print 'Infinity-norm =', x.norm(PETSc.NormType.NORM_INFINITY)
print 'Minimum element in x (index, value) =', x.min()
print 'Maximum element in x (index, value) =', x.max()