Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

mutivec implemented in algorithms, random method modifed in utils #19

Merged
merged 25 commits into from
Apr 10, 2024
Merged

mutivec implemented in algorithms, random method modifed in utils #19

Show file tree
Hide file tree
Changes from 22 commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
1386873
check Borthogonalization operation
Apr 3, 2024
beb97be
parRandom perturb needs to be fixed
Apr 3, 2024
de63f97
extended parRandom to multivector
Apr 3, 2024
5da7f89
multivec added, random modified
Apr 3, 2024
744b849
misc script not needed
Apr 4, 2024
09964ee
misc script removed
Apr 4, 2024
ef31688
changes to init and dot in multivec
Apr 4, 2024
efb1aee
Update random.py
uvilla Apr 4, 2024
34b7cf9
Update random.py
uvilla Apr 4, 2024
cd44b72
minor comment fix to check if CI working
Apr 4, 2024
e27f0a1
attempt_v1 multivector testing
Apr 4, 2024
e17bc7f
attempt_v2 multivector testing
Apr 4, 2024
0656851
partial mypy, multivec, doublePassG
Apr 5, 2024
4ccf3e2
modifications to multivec, Eigensolver
Apr 5, 2024
dfb25cd
eigendecomposition added, need to implement check
Apr 5, 2024
be91584
eigendecompositon implemented, need to modfiy tests and mypy cleanup
Apr 5, 2024
6595a26
test for eigendecomposition added
Apr 5, 2024
9b639e6
CI_check_v1
Apr 5, 2024
6c1175e
CI_check_v2
Apr 5, 2024
7f7af29
delted multivec script from misc
Apr 5, 2024
ace44e0
checking CI post changes v_1
Apr 9, 2024
8ef6223
checking CI post changes v_2
Apr 9, 2024
1902849
corrections to linalg, randomizedeigen solver, multivector testing
Apr 10, 2024
ee5e05e
formatting fix
Apr 10, 2024
9d2273a
remove commented code
Apr 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 13 additions & 2 deletions .github/workflows/CI_testing.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,14 +25,25 @@ jobs:
run: |
pip install ruff && ruff check . && ruff format --check .

- name: multivector testing
run: |
cd ./hippylibX/test &&
mpirun -n 2 python3 test_multivector.py
uvilla marked this conversation as resolved.
Show resolved Hide resolved

- name: eigendecomposition testing
run: |
cd ./hippylibX/test &&
mpirun -n 2 python3 test_eigendecomposition.py

- name: run serial check
run: |
cd ./hippylibX/test &&
mpirun -n 1 python3 testing_suite_file.py
mpirun -n 1 python3 test_model.py

- name: run parallel check
run: |
cd ./hippylibX/test &&
mpirun -n 2 python3 testing_suite_file.py
mpirun -n 2 python3 test_model.py



80 changes: 58 additions & 22 deletions example/poisson_dirichlet_example.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,13 @@
import os
import dolfinx.fem.petsc
from matplotlib import pyplot as plt
from typing import Sequence, Dict


sys.path.append(os.environ.get("HIPPYLIBX_BASE_DIR", "../"))
import hippylibX as hpx


def master_print(comm, *args, **kwargs):
if comm.rank == 0:
print(*args, **kwargs)


class Poisson_Approximation:
def __init__(self, f: float):
self.f = f
Expand All @@ -42,7 +39,9 @@ def __call__(self, u: dlx.fem.Function, m: dlx.fem.Function) -> ufl.form.Form:
return 0.5 / self.sigma2 * ufl.inner(u - self.d, u - self.d) * self.dx


def run_inversion(nx: int, ny: int, noise_variance: float, prior_param: dict) -> None:
def run_inversion(
nx: int, ny: int, noise_variance: float, prior_param: Dict[str, float]
) -> Dict[str, Dict[str, float]]:
sep = "\n" + "#" * 80 + "\n"
comm = MPI.COMM_WORLD
rank = comm.rank
Expand All @@ -56,15 +55,15 @@ def run_inversion(nx: int, ny: int, noise_variance: float, prior_param: dict) ->
Vh_phi.dofmap.index_map.size_global * Vh_phi.dofmap.index_map_bs,
Vh_m.dofmap.index_map.size_global * Vh_m.dofmap.index_map_bs,
]
master_print(comm, sep, "Set up the mesh and finite element spaces", sep)
master_print(comm, "Number of dofs: STATE={0}, PARAMETER={1}".format(*ndofs))
hpx.master_print(comm, sep, "Set up the mesh and finite element spaces", sep)
hpx.master_print(comm, "Number of dofs: STATE={0}, PARAMETER={1}".format(*ndofs))

# dirichlet B.C.
uD = dlx.fem.Function(Vh[hpx.STATE])
uD.interpolate(lambda x: x[1])
uD.x.scatter_forward()

def top_bottom_boundary(x):
def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
return np.logical_or(np.isclose(x[1], 1), np.isclose(x[1], 0))

fdim = msh.topology.dim - 1
Expand Down Expand Up @@ -155,6 +154,17 @@ def top_bottom_boundary(x):

x = solver.solve(x)

if solver.converged:
hpx.master_print(comm, "\nConverged in ", solver.it, " iterations.")
else:
hpx.master_print(comm, "\nNot Converged")

hpx.master_print(
comm, "Termination reason: ", solver.termination_reasons[solver.reason]
)
hpx.master_print(comm, "Final gradient norm: ", solver.final_grad_norm)
hpx.master_print(comm, "Final cost: ", solver.final_cost)

m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")

Expand All @@ -179,16 +189,6 @@ def top_bottom_boundary(x):
) as vtx:
vtx.write(0.0)

if solver.converged:
master_print(comm, "\nConverged in ", solver.it, " iterations.")
else:
master_print(comm, "\nNot Converged")
master_print(
comm, "Termination reason: ", solver.termination_reasons[solver.reason]
)
master_print(comm, "Final gradient norm: ", solver.final_grad_norm)
master_print(comm, "Final cost: ", solver.final_cost)

optimizer_results = {}
if (
solver.termination_reasons[solver.reason]
Expand All @@ -198,10 +198,40 @@ def top_bottom_boundary(x):
else:
optimizer_results["optimizer"] = False

Hmisfit = hpx.ReducedHessian(model, misfit_only=True)

k = 80
p = 20
if rank == 0:
print(
"Double Pass Algorithm. Requested eigenvectors: {0}; Oversampling {1}.".format(
k, p
)
)

temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
Omega = hpx.MultiVector(temp_para_vec, k + p)
temp_para_vec.destroy()

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}

final_results = {
"data_misfit_True": data_misfit_True,
"data_misfit_False": data_misfit_False,
"optimizer_results": optimizer_results,
"eigen_decomposition_results": eigen_decomposition_results,
}

return final_results
Expand All @@ -213,9 +243,15 @@ def top_bottom_boundary(x):
ny = 64
noise_variance = 1e-4
prior_param = {"gamma": 0.03, "delta": 0.3}
run_inversion(nx, ny, noise_variance, prior_param)

final_results = run_inversion(nx, ny, noise_variance, prior_param)
k, d = (
final_results["eigen_decomposition_results"]["k"],
final_results["eigen_decomposition_results"]["d"],
)
comm = MPI.COMM_WORLD
if comm.rank == 0:
plt.savefig("poisson_result_FD_Gradient_Hessian_Check")
plt.show()
plt.figure()
plt.plot(range(0, k), d, "b*", range(0, k), np.ones(k), "-r")
plt.yscale("log")
plt.savefig("poisson_Eigen_Decomposition_results.png")
79 changes: 57 additions & 22 deletions example/poisson_dirichlet_example_reg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,13 @@
import os
import dolfinx.fem.petsc
from matplotlib import pyplot as plt
from typing import Sequence, Dict


sys.path.append(os.environ.get("HIPPYLIBX_BASE_DIR", "../"))
import hippylibX as hpx


def master_print(comm, *args, **kwargs):
if comm.rank == 0:
print(*args, **kwargs)


class Poisson_Approximation:
def __init__(self, f: float):
self.f = f
Expand All @@ -42,7 +39,9 @@ def __call__(self, u: dlx.fem.Function, m: dlx.fem.Function) -> ufl.form.Form:
return 0.5 / self.sigma2 * ufl.inner(u - self.d, u - self.d) * self.dx


def run_inversion(nx: int, ny: int, noise_variance: float, prior_param: dict) -> None:
def run_inversion(
nx: int, ny: int, noise_variance: float, prior_param: Dict[str, float]
) -> Dict[str, Dict[str, float]]:
sep = "\n" + "#" * 80 + "\n"
comm = MPI.COMM_WORLD
rank = comm.rank
Expand All @@ -56,15 +55,15 @@ def run_inversion(nx: int, ny: int, noise_variance: float, prior_param: dict) ->
Vh_phi.dofmap.index_map.size_global * Vh_phi.dofmap.index_map_bs,
Vh_m.dofmap.index_map.size_global * Vh_m.dofmap.index_map_bs,
]
master_print(comm, sep, "Set up the mesh and finite element spaces", sep)
master_print(comm, "Number of dofs: STATE={0}, PARAMETER={1}".format(*ndofs))
hpx.master_print(comm, sep, "Set up the mesh and finite element spaces", sep)
hpx.master_print(comm, "Number of dofs: STATE={0}, PARAMETER={1}".format(*ndofs))

# dirichlet B.C.
uD = dlx.fem.Function(Vh[hpx.STATE])
uD.interpolate(lambda x: x[1])
uD.x.scatter_forward()

def top_bottom_boundary(x):
def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
return np.logical_or(np.isclose(x[1], 1), np.isclose(x[1], 0))

fdim = msh.topology.dim - 1
Expand Down Expand Up @@ -164,6 +163,16 @@ def top_bottom_boundary(x):

x = solver.solve(x)

if solver.converged:
hpx.master_print(comm, "\nConverged in ", solver.it, " iterations.")
else:
hpx.master_print(comm, "\nNot Converged")
hpx.master_print(
comm, "Termination reason: ", solver.termination_reasons[solver.reason]
)
hpx.master_print(comm, "Final gradient norm: ", solver.final_grad_norm)
hpx.master_print(comm, "Final cost: ", solver.final_cost)

m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")

Expand All @@ -190,16 +199,6 @@ def top_bottom_boundary(x):
) as vtx:
vtx.write(0.0)

if solver.converged:
master_print(comm, "\nConverged in ", solver.it, " iterations.")
else:
master_print(comm, "\nNot Converged")
master_print(
comm, "Termination reason: ", solver.termination_reasons[solver.reason]
)
master_print(comm, "Final gradient norm: ", solver.final_grad_norm)
master_print(comm, "Final cost: ", solver.final_cost)

optimizer_results = {}
if (
solver.termination_reasons[solver.reason]
Expand All @@ -209,10 +208,40 @@ def top_bottom_boundary(x):
else:
optimizer_results["optimizer"] = False

Hmisfit = hpx.ReducedHessian(model, misfit_only=True)

k = 80
p = 20
if rank == 0:
print(
"Double Pass Algorithm. Requested eigenvectors: {0}; Oversampling {1}.".format(
k, p
)
)

temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
Omega = hpx.MultiVector(temp_para_vec, k + p)
temp_para_vec.destroy()

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}

final_results = {
"data_misfit_True": data_misfit_True,
"data_misfit_False": data_misfit_False,
"optimizer_results": optimizer_results,
"eigen_decomposition_results": eigen_decomposition_results,
}

return final_results
Expand All @@ -225,9 +254,15 @@ def top_bottom_boundary(x):
ny = 64
noise_variance = 1e-4
prior_param = {"gamma": 0.08, "delta": 0.8}
run_inversion(nx, ny, noise_variance, prior_param)

final_results = run_inversion(nx, ny, noise_variance, prior_param)
k, d = (
final_results["eigen_decomposition_results"]["k"],
final_results["eigen_decomposition_results"]["d"],
)
comm = MPI.COMM_WORLD
if comm.rank == 0:
plt.savefig("poisson_result_FD_Gradient_Hessian_Check")
plt.show()
plt.figure()
plt.plot(range(0, k), d, "b*", range(0, k), np.ones(k), "-r")
plt.yscale("log")
plt.savefig("poisson_Eigen_Decomposition_results.png")
Loading
Loading