From e62698e1bbc3c6493ad1bcfb3971ba86a244adc7 Mon Sep 17 00:00:00 2001 From: David Seus <david.seus@ians.uni-stuttgart.de> Date: Wed, 3 Apr 2019 10:21:37 +0200 Subject: [PATCH] rewrite self._eval_sources() to operate on a single domain --- LDDsimulation/LDDsimulation.py | 94 +++++++++++++++++++--------------- 1 file changed, 54 insertions(+), 40 deletions(-) diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py index f6ed481..3da0268 100644 --- a/LDDsimulation/LDDsimulation.py +++ b/LDDsimulation/LDDsimulation.py @@ -145,9 +145,6 @@ class LDDsimulation(object): self._init_interfaces() self._init_subdomains() self._init_initial_values() - self._eval_sources() - # init Dirichlet Boundary values for the first time - self.set_DirichletBC(time = 0, first_init = True) def LDDsolver(self, time: float): """ calulate the solution on all patches for the next timestep @@ -165,6 +162,11 @@ class LDDsimulation(object): # and for ind, subdomain in self.subdomain.items(): subdomain_calc_isdone.update({ind, False}) + # update the time of the source terms before starting the iteration. + # The sources and sinks are the same throughout the iteration. + self._eval_sources(interpolation_degree = 2, + time = time, + subdomain_index = ind) # the following only needs to be done when time is not equal 0 since # our initialisation functions already took care of setting what follows # for the initial iteration step. @@ -177,6 +179,7 @@ class LDDsimulation(object): # TODO recheck if subdomain.pressure_prev_iter = subdomain.pressure # needs to be set also + ### actual iteration starts here # gobal stopping criterion for the iteration. all_subdomains_done = False @@ -189,12 +192,12 @@ class LDDsimulation(object): # finished. if not subdomain_calc_isdone[sd_index]: # solve the problem on subdomain - self.Lsolver_step(subdomain = subdomain,# + self.Lsolver_step(subdomain_index = sd_index,# iteration = iteration,# debug = True ) subsequent_iterations_err = self.calc_iteration_error( - subdomain = subdomain,# + subdomain_index = sd_index,# error_type = "L2" ) # check if error criterion has been met @@ -226,6 +229,26 @@ class LDDsimulation(object): iteration += 1 # end iteration while loop. + def Lsolver_step(subdomain_index: int],# + iteration: int,# + debug: bool = False) -> None: + """ L-scheme solver iteration step for an object of class subdomain + + Lsolver_step implements L-scheme solver iteration step for an subdomain + with index subdomain_index, i.e. for the model form (already in L-scheme + form) stored in the domain patch subdomain. + + # parameters + subdomain_index: int subdomain index that defines + which subdomain to perform + the iteration on. + iteration: int iteration number + debug: bool flag to toggle weather or not + to write out each iteration. + """ + + + ## Private methods def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,# mesh_resolution: float = None) -> None: @@ -444,43 +467,34 @@ class LDDsimulation(object): {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},# ) - def _eval_sources(self, interpolation_degree: int = 2, time: float = 0): + def _eval_sources(self, interpolation_degree: int = 2, # + time: float = 0, + subdomain_index: int): """ evaluate time dependent source terms or initialise them if time == 0 """ + subdomain = self.subdomain[subdomain_index] if np.abs(time) < self.tol: # here t = 0 and we have to initialise the sources. - for num, subdomain in self.subdomain.items(): - mesh = subdomain.mesh - V = subdomain.function_space - # p0 contains both pw_0 and pnw_0 for subdomain[num] - f = self.sources[num] - if subdomain.isRichards: - # note that the default interpolation degree is 2 - fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time) - # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])} - subdomain.source = {'wetting' : fw} - else: - fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time) - fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time) - # subdomain.pressure_prev_iter = {'wetting' : df.interpolate(fw, V['wetting']),# - # 'nonwetting' : df.interpolate(fnw, V['nonwetting'])} - # subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(fw, V['wetting']),# - # 'nonwetting' : df.interpolate(fnw, V['nonwetting'])} - subdomain.source = {'wetting' : fw,# - 'nonwetting' : fnw} + mesh = subdomain.mesh + V = subdomain.function_space + # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index] + f = self.sources[subdomain_index] + if subdomain.isRichards: + # note that the default interpolation degree is 2 + fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time) + # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])} + subdomain.source = {'wetting' : fw} + else: + fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time) + fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time) + subdomain.source = {'wetting' : fw,# + 'nonwetting' : fnw} else: - for num, subdomain in self.subdomain.items(): - # mesh = subdomain.mesh - # V = subdomain.function_space - # # p0 contains both pw_0 and pnw_0 for subdomain[num] - # f = self.sources[num] - if subdomain.isRichards: - # note that the default interpolation degree is 2 - # fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time) - # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])} - subdomain.source['wetting'].t = time - else: - # fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time) - # fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time) - subdomain.source['wetting'].t = time - subdomain.source['nonwetting'].t = time + # here the sources are already initiated and stored in the subdomain. + # they must be updated in time. + if subdomain.isRichards: + # note that the default interpolation degree is 2 + subdomain.source['wetting'].t = time + else: + subdomain.source['wetting'].t = time + subdomain.source['nonwetting'].t = time -- GitLab