Skip to content

Commit

Permalink
MAINT: Switch to relative criteria for Bregman and Newton.
Browse files Browse the repository at this point in the history
  • Loading branch information
@jwboth
jwboth committed Sep 9, 2024
1 parent 3626869 commit eec7483
Showing 1 changed file with 29 additions and 26 deletions.
55 changes: 29 additions & 26 deletions src/darsia/measure/wasserstein.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1556,7 +1556,7 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
# - residual
if self.verbose:
print(
"Newton iter. \t| W^1 \t\t| Δ W^1 \t| Δ flux \t| residual",
"Newton iter. \t| W^1 \t\t| Δ W^1 / W^1 \t| Δ flux \t| residual",
"\n",
"""---------------|---------------|---------------|---------------|"""
"""---------------""",
Expand All @@ -1569,8 +1569,8 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
try:
# Keep track of old flux, and old distance
old_solution_i = solution_i.copy()
old_flux = solution_i[self.flux_slice]
old_distance = self.l1_dissipation(old_flux)
flux = solution_i[self.flux_slice]
old_distance = self.l1_dissipation(flux)

# Assemble linear problem in Newton step
tic = time.time()
Expand Down Expand Up @@ -1603,8 +1603,8 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
stats_i["time_acceleration"] = time.time() - tic

# Update discrete W1 distance
new_flux = solution_i[self.flux_slice]
new_distance = self.l1_dissipation(new_flux)
flux = solution_i[self.flux_slice]
new_distance = self.l1_dissipation(flux)

# Update increment
increment = solution_i - old_solution_i
Expand Down Expand Up @@ -1637,7 +1637,9 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
# - flux increment
# - residual
if self.verbose:
distance_increment = convergence_history["distance_increment"][-1]
distance_increment = (
convergence_history["distance_increment"][-1] / new_distance
)
flux_increment = (
convergence_history["flux_increment"][-1]
/ convergence_history["flux_increment"][0]
Expand Down Expand Up @@ -1827,7 +1829,7 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
# Print header
if self.verbose:
print(
"Bregman iter. \t| W^1 \t\t| Δ W^1 \t| Δ aux/force \t| mass residual",
"Bregman iter. \t| W^1 \t\t| Δ W^1/W^1 \t| Δ aux/force \t| mass residual",
"\n",
"""---------------|---------------|---------------|---------------|"""
"""---------------""",
Expand All @@ -1848,12 +1850,11 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
)

# Initialize Bregman variables
old_flux = solution_i[self.flux_slice]
old_aux_flux = self._shrink(old_flux, shrink_factor)
old_force = old_flux - old_aux_flux
old_distance = self.l1_dissipation(old_flux)
flux = solution_i[self.flux_slice]
old_aux_flux = self._shrink(flux, shrink_factor)
old_force = flux - old_aux_flux
old_distance = self.l1_dissipation(flux)

new_flux = old_flux.copy()
iter = 0

# Control the update of the Bregman weight
Expand All @@ -1874,7 +1875,7 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
l_scheme_mixed_darcy,
weight,
shrink_factor,
) = self._update_regularization(old_flux, bregman_homogeneous)
) = self._update_regularization(flux, bregman_homogeneous)
# 1. Make relaxation step (solve quadratic optimization problem)
# Here, re-initialize the aux flux and force with zero values again.
rhs_i = rhs.copy()
Expand All @@ -1886,17 +1887,18 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
rhs_i,
reuse_solver=False,
)
new_flux = solution_i[self.flux_slice]
flux = solution_i[self.flux_slice]
stats_i["time_solve"] = time.time() - tic
stats_i["time_assemble"] = time_assemble

# 2. Shrink step for vectorial fluxes.
tic = time.time()
new_aux_flux = self._shrink(new_flux, shrink_factor)
new_aux_flux = self._shrink(flux, shrink_factor)
stats_i["time_shrink"] = time.time() - tic

# 3. Update force
new_force = new_flux - new_aux_flux
tic = time.time()
new_force = flux - new_aux_flux

else:
# 1. Make relaxation step (solve quadratic optimization problem)
Expand All @@ -1913,17 +1915,18 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
rhs_i,
reuse_solver=iter > 0,
)
new_flux = solution_i[self.flux_slice]
flux = solution_i[self.flux_slice]
stats_i["time_solve"] = time.time() - tic
stats_i["time_assemble"] = time_assemble

# 2. Shrink step for vectorial fluxes.
tic = time.time()
new_aux_flux = self._shrink(new_flux + old_force, shrink_factor)
new_aux_flux = self._shrink(flux + old_force, shrink_factor)
stats_i["time_shrink"] = time.time() - tic

# 3. Update force
new_force = old_force + new_flux - new_aux_flux
tic = time.time()
new_force = old_force + flux - new_aux_flux

# Apply Anderson acceleration to flux contribution (the only nonlinear part).
tic = time.time()
Expand All @@ -1938,15 +1941,15 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
stats_i["time_acceleration"] = time.time() - tic

# Update distance
new_distance = self.l1_dissipation(new_flux)
new_distance = self.l1_dissipation(flux)

# Determine the error in the mass conservation equation
mass_conservation_residual = (
self.div.dot(new_flux) - rhs[self.pressure_slice]
self.div.dot(flux) - rhs[self.pressure_slice]
)

# Reference values
flux_ref = np.linalg.norm(new_flux, 2)
flux_ref = np.linalg.norm(flux, 2)
mass_ref = np.linalg.norm(rhs[self.pressure_slice], 2)

# Determine increments
Expand Down Expand Up @@ -1986,7 +1989,7 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
)
aux_force_increment = (
convergence_history["aux_force_increment"][-1]
# / convergence_history["aux_force_increment"][0]
/ convergence_history["aux_force_increment"][0]
)
mass_conservation_residual = convergence_history[
"mass_conservation_residual"
Expand All @@ -2013,6 +2016,7 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
< tol_increment
* convergence_history["aux_force_increment"][0]
and convergence_history["distance_increment"][-1]
/ new_distance
< tol_distance
and convergence_history["mass_conservation_residual"][-1]
< tol_residual
Expand All @@ -2021,7 +2025,6 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
break

# Update Bregman variables
old_flux = new_flux.copy()
old_aux_flux = new_aux_flux.copy()
old_force = new_force.copy()
old_distance = new_distance
Expand All @@ -2031,9 +2034,9 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
break

# Solve for the pressure by solving a single Newton iteration
newton_jacobian, _, _ = self._update_regularization(new_flux)
newton_jacobian, _, _ = self._update_regularization(flux)
solution_i = np.zeros_like(rhs)
solution_i[self.flux_slice] = new_flux.copy()
solution_i[self.flux_slice] = flux.copy()
newton_residual = self.optimality_conditions(rhs, solution_i)
newton_update, _ = self.linear_solve(
newton_jacobian, newton_residual, solution_i
Expand Down

0 comments on commit eec7483

Please sign in to comment.