diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 79d657b237bc504f49f9a8e2e820781472fcf0b5..a31b04fc4825f4f803fa7df8a04edb5da7fb80f2 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -468,6 +468,7 @@ class LDDsimulation(object):
- the calculated solution pressure from the previous timestep
is assigned to pressure_prev for each subdomain.
- gl0 terms are calculated.
+ - write pressure to previous iteration dictionaries on all interfaces.
PARAMETERS
------------------------------------
@@ -532,7 +533,12 @@ class LDDsimulation(object):
if debug:
print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure[{phase}])", "red"), " =\n", f"{id(subdomain.pressure[phase])}")
print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
-
+ # we need to update the previous pressure dictionaries of all interfaces
+ # of the subdomain with the current solution, so that the subdomain.calc_gli_term
+ # method works properly. The current pressure dictionaries should already
+ # be populated with the current solution, sinces the solver writes the pressure
+ # to interfaces after each iteration.
+ subdomain.write_pressure_to_interfaces(previous_iter=True)
subdomain.calc_gl0_term()
return solution_over_iteration_within_timestep
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 76832d4290b1192f4cc1f3021ae8492c792927c1..564f415468aa0a8ca49677fae67afa41b936fc3a 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -471,6 +471,7 @@ class DomainPatch(df.SubDomain):
Caution: The array self.gli_assembled needs to be added to the rhs with
a minus sign by the solver!
"""
+ # this is the number of the current iteration of the subdomain.
iteration = self.iteration_number
subdomain = self.subdomain_index
dt = self.timestep_size
@@ -489,6 +490,10 @@ class DomainPatch(df.SubDomain):
interface = self.interface[ind]
# neighbour index
neighbour = interface.neighbour[subdomain]
+ # update the current_iteration number of subdomain stored in the interface
+ # to the current iteration number of the subdomain.
+ interface.current_iteration[subdomain] = iteration
+ # save the iteration number of the neighbour
neighbour_iter_num = interface.current_iteration[neighbour]
ds = self.ds(interface.marker_value)
# needed for the read_gli_dofs() functions
@@ -527,37 +532,26 @@ class DomainPatch(df.SubDomain):
# read previous neighbouring pressure values from interface and write them to
# the function pa_prev. We need to make a distinction first.
- if neighbour_iter_num <= iteration:
- # in case neighbour_iter_num == iteration the neighbour
- # has not yet iterated beyond
+ if iteration == neighbour_iter_num + 1:
+ # in this case the neighbour has not yet iterated beyond
# the iteration we are currently in, so
- # interface.pressure_values holds the "previous pressure"
- # we need to read from the interface to calculate gli.
- # In case neighbour_iter_num < iteration the neighbour has
- # already reached the error tolerance and has stopped iterating.
- # We therefor read again the current and not the previous pressure
+ # interface.pressure_values[neighbour] holds the "previous pressure"
+ # we need to read from the interface to calculate the current gli.
take_previous_pressure_from_interface = False
- if neighbour_iter_num + 1 < iteration:
- print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
- f"and iteration = {iteration} on subdomain {subdomain}\n",
- "It is expected that this case should only be happening if the \n",
- f"neighbour (index {neighbour}) has already finished its iteration, i.e. reached\n",
- f"the error criterion. In this case we should have {iteration} > {neighbour_iter_num}.\n",
- "If this is not the case revisite DomainPatch.calc_gli_term().")
- else:
- # the case we think here is neighbour_iter_num = iteration + 1.
- # in this case the neighbour has iterated once beyond
- # the iteration we are currently in, so
- # interface.pressure_values_prev holds the "previous pressure"
+
+ elif iteration == neighbour_iter_num:
+ # this is the only possible other case. In this case the
+ # neighbour has iterated once beyond the iteration we are
+ # currently in, so interface.pressure_values_prev holds the "previous pressure"
# we need to read from the interface to calculate gli.
- # neighbour_iter_num > iteration + 1 should not appear since
- # we only stop to iterate if we reach error tolerance and no
- # subdomain is left out while iterating.
take_previous_pressure_from_interface = True
- if neighbour_iter_num > iteration + 1:
- print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
- f"and iteration = {iteration} on subdomain {subdomain}\n",
- "This case should not be happening.")
+ else:
+ raise RuntimeError("\n!!!!!!!!! Warning !!!!!!!!!!!!!! \n"+\
+ "this case should not appear\n"+"neighbour_iter_num =" + "{}".format(neighbour_iter_num)\
+ + "for neighbour = "+ "{}".format(neighbour)+\
+ "and iteration = " + "{} on subdomain".format(iteration)+ "{}".format(subdomain)\
+ + "You suck! Revisite DomainPatch.calc_gli_term() and LDDsimulation.LDDsolver ASAP.")
+
pa_prev = df.interpolate(df.Constant(0), self.function_space[phase])
interface.write_pressure_dofs(to_function = pa_prev,#self.neighbouring_interface_pressure[phase], #
interface_dofs = interface_dofs[phase],#
@@ -572,8 +566,8 @@ class DomainPatch(df.SubDomain):
# pa_prev = self.neighbouring_interface_pressure[phase]
phi_a = self.testfunction[phase]
gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
- gli_p_prev_part = gli_pressure_part.get_local()
- gli = gli_p_prev_part - gli_prev_of_neighbour
+ gli_pressure_part_array = gli_pressure_part.get_local()
+ gli = gli_pressure_part_array - gli_prev_of_neighbour
# check if there are shared dofs with another interface.
if dofs_in_common_with_other_interfaces[phase].size != 0:
for common_dof in dofs_in_common_with_other_interfaces[phase]:
@@ -599,7 +593,7 @@ class DomainPatch(df.SubDomain):
print(f"before writing to interface[{ind}] gli_assembled[{phase}]=\n",self.gli_assembled[phase])
# write gli dofs to interface dicts for communiction. We have to
# distinguish different cases.
- if neighbour_iter_num == iteration:
+ if iteration == neighbour_iter_num + 1:
# in this case, the neighbour has not done the next iteration
# and will need both pressure and gli_term_prev of the current
# subdomain. Therefor, we save the the calculated gli_terms to
@@ -613,7 +607,8 @@ class DomainPatch(df.SubDomain):
subdomain_ind = subdomain,#
previous_iter = False
)
- if neighbour_iter_num -1 == iteration:
+ else:
+ # This should be the case iteration == neighbour_iter_num:
# in this case, the neighbour has done the next iteration already
# and will not need both pressure_prev and gli_term_prev of the current
# subdomain. Therefor, we save the the calculated gli_terms to
@@ -638,7 +633,6 @@ class DomainPatch(df.SubDomain):
interface.pressure_values_prev[neighbour][phase].update(
{phase: interface.pressure_values[neighbour][phase].copy()}
)
- interface.current_iteration[subdomain] += 1
### END calc_gli_term
def null_all_interface_iteration_numbers(self) -> None: