fix PCP

pylbm_ui/responses.py

@@ -12,6 +12,7 @@
 import copy
 import pylbm
 
+from .widgets.debug import debug
 class FromConfig:
     pass
 
@@ -23,7 +24,7 @@ class AfterSimulation:
     pass
 
 class Sigma(FromConfig):
-    def __init__(self, s, log10=True):
+    def __init__(self, s, log10=False):
         self.log10 = log10
         self.s = s
 
@@ -36,8 +37,7 @@ def __call__(self, config, extra_config=None):
         return np.log10(output) if self.log10 else output
 
 class S(FromConfig):
-    def __init__(self, s, log10=True):
-        self.log10 = log10
+    def __init__(self, s):
         self.s = s
 
     def __call__(self, config, extra_config=None):
@@ -46,10 +46,10 @@ def __call__(self, config, extra_config=None):
         params.update(extra_config)
 
         output = params[self.s]
-        return np.log10(output) if self.log10 else output
+        return  output
 
 class Diff(FromConfig):
-    def __init__(self, s, with_dx=True, log10=True):
+    def __init__(self, s, with_dx=True, log10=False):
         self.log10 = log10
         self.with_dx = with_dx
         self.s = s
@@ -89,7 +89,7 @@ def __call__(self, config, extra_config=None):
         return True
 
 class Error(AfterSimulation):
-    def __init__(self, ref_solution, expr, log10=False, relative=False):
+    def __init__(self, ref_solution, expr, log10=True, relative=False):
         self.ref_solution = ref_solution
         self.expr = expr
         self.log10 = log10
@@ -117,46 +117,60 @@ def __call__(self, sol):
             norm /= np.linalg.norm(self.ref_solution)
         return np.log10(norm) if self.log10 else norm
 
+@debug
 class ErrorStd(DuringSimulation):
-    def __init__(self, field, ref_func, expr, call_at=10, log10=True):
+    def __init__(self, field, ref_func, expr, call_at=0.92, log10=True):
         self.field = field
         self.ref_func = ref_func
         self.expr = expr
         self.log10 = log10
         self.error = []
         self.nite = 0
         self.call_at = call_at
+        #print('hello')
 
-    def __call__(self, sol):
+    def __call__(self, sol, duration):
         self.nite += 1
-        if self.nite == self.call_at:
+        startTime = duration*self.call_at
+        startIt = int(startTime/sol.dt)
+
+        if sol.t >= startTime:
+        #if self.nite == startIt:
+            #print(sol.t, duration, startTime)
             self.nite = 0
             domain = sol.domain
             time_e = sol.t
             ref_solution = self.ref_func(time_e, domain.x, field=self.field)
 
-            func = sp.lambdify(list(self.expr.atoms(sp.Symbol)), self.expr, "numpy", dummify=False)
-            to_subs = {str(k): sol.m[k] for k in sol.scheme.consm.keys()}
-            to_subs.update({str(k): v for k, v in sol.scheme.param.items()})
+            #func = sp.lambdify(list(self.expr.atoms(sp.Symbol)), self.expr, "numpy", dummify=False)
+            #to_subs = {str(k): sol.m[k] for k in sol.scheme.consm.keys()}
+            #to_subs.update({str(k): v for k, v in sol.scheme.param.items()})
+
+            #args = {str(s): to_subs[str(s)] for s in self.expr.atoms(sp.Symbol)}
+            #data = func(**args)
+            ## remove element with NaN values
+            #solid_cells = sol.domain.in_or_out != sol.domain.valin
 
-            args = {str(s): to_subs[str(s)] for s in self.expr.atoms(sp.Symbol)}
-            data = func(**args)
-            # remove element with NaN values
-            solid_cells = sol.domain.in_or_out != sol.domain.valin
+            #vmax = sol.domain.stencil.vmax
+            #ind = []
+            #for vm in vmax:
+                #ind.append(slice(vm, -vm))
+            #ind = np.asarray(ind)
 
-            vmax = sol.domain.stencil.vmax
-            ind = []
-            for vm in vmax:
-                ind.append(slice(vm, -vm))
-            ind = np.asarray(ind)
+            #data[solid_cells[tuple(ind)]] = 0
+            #norm = np.linalg.norm(ref_solution - data)
+
+            error = Error(ref_solution, self.expr, log10=False, relative=False)
+            norm = error(sol)
 
-            data[solid_cells[tuple(ind)]] = 0
-            norm = np.linalg.norm(ref_solution - data)
             self.error.append(norm)
+            #print (len(self.error))
 
     def value(self):
         std = np.std(np.asarray(self.error))
         return np.log10(std) if self.log10 else std
+        #std = len(self.error)
+        #return std
 
 class ErrorAvg(ErrorStd):
     def value(self):
@@ -236,12 +250,4 @@ def __call__(self, sol):
         fig.savefig(self.filename, dpi=300)
 
 
-class StdError(AfterSimulation):
-    def __init__(self, ref_solution, field):
-        self.ref_solution = ref_solution
-        self.field = field
-
-    def __call__(self, simulation):
-        data = simulation.get_data(self.field, 0)
-        return np.linalg.norm(self.ref_solution - data)
 

pylbm_ui/widgets/parametric_study.py

@@ -46,39 +46,43 @@ def run_simulation(args):
         if isinstance(r, FromConfig):
             output[i] = r(simu_cfg, sample)
 
-    def test_nan():
-        c = list(sol.scheme.consm.keys())[0]
+    def test_unstab():
+        c = list(sol.scheme.consm.keys())[0] ##??? are we sure that it is the mass field?
         sol.m_halo[c][solid_cells] = 0
         data = sol.m[c]
-        if np.isnan(np.sum(data)) or np.any(np.abs(data)>1e20):
+        if np.isnan(np.sum(data)) or np.any(np.abs(data)>1e10):
+            return True
+        if np.any(data<=0.): ## avoid negative mass
             return True
         return False
+
 
     actions = [r for r in responses if isinstance(r, DuringSimulation)]
 
-    nan_detected = False
+    unstable = False
     nite = 0
-    while sol.t <= duration and not nan_detected:
+    niteStab =  int(duration/sol.dt/10) # the number of stability check during the simulation = 10 
+    while sol.t <= duration and not unstable:
         sol.one_time_step()
 
         for a in actions:
-            a(sol)
+            a(sol, duration)
 
-        if nite == 200:
+        if nite == niteStab:
             nite = 0
-            nan_detected = test_nan()
+            unstable = test_unstab()
         nite += 1
 
-    nan_detected |= test_nan()
-
-    for i, r in enumerate(responses):
-        if isinstance(r, AfterSimulation):
-            output[i] = r(sol)
-        elif isinstance(r, DuringSimulation):
-            output[i] = r.value()
+    unstable |= test_unstab()
 
+    if not unstable: # avoid meaningless responses values (or empty plots) when simulation is unstable 
+        for i, r in enumerate(responses):
+            if isinstance(r, AfterSimulation):
+                output[i] = r(sol)
+            elif isinstance(r, DuringSimulation):
+                output[i] = r.value()
 
-    return [not nan_detected] + output
+    return [not unstable] + output
 
 skopt_method = {'Latin hypercube': Lhs,
                 'Sobol': Sobol,
@@ -172,15 +176,15 @@ def start_PS(self, widget, event, data):
             self.dialog.check_path(path)
 
             asyncio.ensure_future(self.run_study(path))
-            # self.run_study(path)
+            self.run_study(path)
         else:
             self.stop_simulation(None)
 
     async def run_study(self, path):
-    # def run_study(self, path):
+    #def run_study(self, path):
         """
-        Create the sampling using the design space defined in the menu and start the
-        parametric study on each sample.
+        Create the sampling using the design space defined in the menu 
+        and run the evaluation of each the each sample point in parallel.
         Then the results is represented by plotly.
         """
         while self.dialog.v_model:
@@ -258,7 +262,9 @@ async def run_study(self, path):
 
             def run_parametric_study():
                 from pathos.multiprocessing import ProcessingPool
-                pool = pp.ProcessPool()
+                import multiprocessing
+                nNodes = int(multiprocessing.cpu_count()/2)
+                pool = pp.ProcessPool(nodes=nNodes)
                 output = pool.map(run_simulation, args)
 
                 dimensions = [dict(values=np.asarray([o[0] for o in output], dtype=np.float64), label='stability')]

pylbm_ui/widgets/responses.py

@@ -25,7 +25,7 @@ def CFL(test_case):
     return pylbm_responses.CFL(eq.rho, vel)
 
 class Error:
-    def __init__(self, field, expr, relative=False, log10=False):
+    def __init__(self, field, expr, relative=False, log10=True):
         self.field = field
         self.expr = expr
         self.relative = relative
@@ -72,16 +72,19 @@ def update_responses(change):
             for name, expr in fields.items():
                 self.responses[f'plot {name}'] = Plot(name, expr)
                 if hasattr(test_case, 'ref_solution'):
-                    self.responses[f'error on {name}'] = Error(name, expr)
-                    self.responses[f'error avg on {name}'] = pylbm_responses.ErrorAvg(name, test_case.ref_solution, expr)
-                    self.responses[f'error std on {name}'] = pylbm_responses.ErrorStd(name, test_case.ref_solution, expr)
-                    self.responses[f'relative error on {name}'] = Error(name, expr, relative=True)
+                    self.responses[f'err_{name}'] = Error(name, expr)
+                    self.responses[f'errAvg_{name}'] = pylbm_responses.ErrorAvg(name, test_case.ref_solution, expr)
+                    self.responses[f'errStd_{name}'] = pylbm_responses.ErrorStd(name, test_case.ref_solution, expr)
+                    self.responses[f'errRel_{name}'] = Error(name, expr, relative=True)
 
             def add_relax(v):
                 self.responses[k] = pylbm_responses.S(v.symb)
-                self.responses[f'sigma for {k}'] = pylbm_responses.Sigma(v.symb)
-                self.responses[f'diff for {k}'] = pylbm_responses.Diff(v.symb)
-                self.responses[f'diff with dx=1 for {k}'] = pylbm_responses.Diff(v.symb, with_dx=False)
+                self.responses[f'sigma_{k[2:]}'] = pylbm_responses.Sigma(v.symb)
+                self.responses[f'diff_{k[2:]}'] = pylbm_responses.Diff(v.symb)
+                self.responses[f'diffOdx_{k[2:]}'] = pylbm_responses.Diff(v.symb, with_dx=False)
+                self.responses[f'LogSigma_{k[2:]}'] = pylbm_responses.Sigma(v.symb, log10=True)
+                self.responses[f'LogDiff_{k[2:]}'] = pylbm_responses.Diff(v.symb, log10=True)
+                self.responses[f'LogDiffOdx_{k[2:]}'] = pylbm_responses.Diff(v.symb, with_dx=False, log10=True)
 
             for k, v in lb_case.__dict__.items():
                 if isinstance(v, RelaxationParameterFinal):
@@ -110,4 +113,4 @@ def get_list(self, path, test_case, simu_cfg):
                 output.append(self.responses[v])
             else:
                 output.append(self.responses[v](path, test_case, simu_cfg))
-        return output
+        return output

schema/Dimension1/Euler/D1q333.py

@@ -65,9 +65,9 @@ def get_dictionary(self):
         )
         # add numerical diffusion
         addvisc = self.addvisc
-        w0eq = (1-addvisc)*w0eq + addvisc/2*la**2*rho
-        w1eq = (1-addvisc)*w1eq + addvisc/2*la**2*q
-        w2eq = (1-addvisc)*w2eq + addvisc/2*la**2*E
+        w0eq = (1-addvisc)*w0eq + addvisc/2*la_**2*rho
+        w1eq = (1-addvisc)*w1eq + addvisc/2*la_**2*q
+        w2eq = (1-addvisc)*w2eq + addvisc/2*la_**2*E
 
         X = sp.symbols('X')
 
@@ -153,7 +153,7 @@ def get_boundary(self):
         }
 
 class D1Q333(D1Q333_general):
-    addvisc = 0.5
+    addvisc = 0.25
     name = 'D1Q333_0'
     tex_name = r'$D_1Q_{{333}}0$'