Allow n instead of 1 4el-WK models; add option to ramp a DBC profile …

…from file by time; bump version

pyproject.toml

@@ -3,7 +3,7 @@ requires = ["setuptools>=61.0.0", "wheel"]
 
 [project]
 name = "ambit-fe"
-version = "v1.2.6"
+version = "v1.2.7"
 description = "A FEniCS-based cardiovascular multi-physics solver"
 authors = [{name = "Marc Hirschvogel", email = "[email protected]"}]
 license = {file = "LICENSE"}

src/ambit_fe/ale/ale_kinematics_constitutive.py

@@ -60,6 +60,10 @@ def stress(self, d_):
 
                 stress += mat.diffusion(matparams_m)
 
+            elif matlaw == 'diffusion_sym':
+
+                stress += mat.diffusion_sym(matparams_m)
+
             else:
 
                 raise NameError('Unknown ALE material law!')

src/ambit_fe/ale/ale_material.py

@@ -69,3 +69,10 @@ def diffusion(self, params):
         D = params['D']
 
         return D*ufl.grad(self.d)
+
+
+    def diffusion_sym(self, params):
+
+        D = params['D']
+
+        return D*ufl.sym(ufl.grad(self.d))

src/ambit_fe/boundaryconditions.py

@@ -80,7 +80,27 @@ def dirichlet_bcs(self, bcdict, V):
                 fle = d['file'] # a single file
                 try: ftype = d['ftype']
                 except: ftype = 'id_val'
-                self.io.readfunction(func, fle, filetype=ftype)
+                # to ramp the file by a time curve
+                try: ramp_curve = d['ramp_curve']
+                except: ramp_curve = None
+                if ramp_curve is not None:
+                    func_ramp, func_file = fem.Function(V), fem.Function(V)
+                    # first read file into function
+                    self.io.readfunction(func_file, fle, filetype=ftype)
+                    # now store ramp curve into function
+                    load_ = expression.template_vector(dim=self.dim)
+                    load_.val_x, load_.val_y, load_.val_z = self.ti.timecurves(d['ramp_curve'])(self.ti.t_init), self.ti.timecurves(d['ramp_curve'])(self.ti.t_init), self.ti.timecurves(d['ramp_curve'])(self.ti.t_init)
+                    func_ramp.interpolate(load_.evaluate)
+                    func_ramp.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
+                    self.ti.funcs_to_update_vec.append({func_ramp : [self.ti.timecurves(d['ramp_curve']), self.ti.timecurves(d['ramp_curve']), self.ti.timecurves(d['ramp_curve'])],
+                                                        'funcs_mult' : [func_file, func]})
+                    self.ti.funcs_to_update_vec_old.append({None : -1}) # DBCs don't need an old state
+                    # now multiply
+                    func.vector.pointwiseMult(func_ramp.vector, func_file.vector)
+                    func.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
+                else:
+                    # read file into function
+                    self.io.readfunction(func, fle, filetype=ftype)
             else:
                 raise RuntimeError("Need to have 'curve', 'val', 'expression', or 'file' specified!")
 

src/ambit_fe/flow0d/cardiovascular0D_2elwindkessel.py

@@ -71,7 +71,6 @@ def setup_arrays(self):
 
             self.v_ids.append(n)
             self.c_ids.append(n)
-
             if self.cq[n] == 'volume':
                 assert(self.vq[n]=='pressure')
                 self.switch_V.append(1), self.cname.append('V'), self.vname.append('p')
@@ -81,7 +80,7 @@ def setup_arrays(self):
             elif self.cq[n] == 'pressure':
                 if self.vq[n] == 'flux':
                     self.switch_V.append(0), self.cname.append('p'), self.vname.append('Q')
-                elif self.vq[0] == 'volume':
+                elif self.vq[n] == 'volume':
                     self.switch_V.append(1), self.cname.append('p'), self.vname.append('V')
                 else:
                     raise ValueError("Unknown variable quantity!")

src/ambit_fe/flow0d/cardiovascular0D_4elwindkesselLpZ.py

@@ -31,14 +31,29 @@ def __init__(self, params, cq, vq, init=True, ode_par=False, comm=None):
         # initialize base class
         super().__init__(init=init, ode_par=ode_par, comm=comm)
 
-        self.num_models = 1 # up to now hard-set to 1
+        # number of (independent) models
+        try: self.num_models = params['num_models']
+        except: self.num_models = 1
+
+        self.R, self.C, self.Z, self.L, self.p_ref = [], [], [], [], []
 
         # parameters
-        self.R = params['R']
-        self.C = params['C']
-        self.Z = params['Z']
-        self.L = params['L']
-        self.p_ref = params['p_ref']
+        for n in range(self.num_models):
+            # resistance
+            try: self.R.append(params['R'+str(n+1)])
+            except: self.R.append(params['R'])
+            # compliance
+            try: self.C.append(params['C'+str(n+1)])
+            except: self.C.append(params['C'])
+            # impedance
+            try: self.Z.append(params['Z'+str(n+1)])
+            except: self.Z.append(params['Z'])
+            # inertance
+            try: self.L.append(params['L'+str(n+1)])
+            except: self.L.append(params['L'])
+            # downstream reference pressure
+            try: self.p_ref.append(params['p_ref'+str(n+1)])
+            except: self.p_ref.append(params['p_ref'])
 
         self.cq = cq
         self.vq = vq
@@ -61,91 +76,117 @@ def setup_arrays(self):
         # number of degrees of freedom - 4 per model
         self.numdof = 4*self.num_models
 
-        self.v_ids = [0]
-
-        if self.cq[0] == 'volume':
-            self.c_ids = [2]
-            assert(self.vq[0]=='pressure')
-            self.switch_V, self.cname, self.vname = 1, 'V', 'p'
-        elif self.cq[0] == 'flux':
-            self.c_ids = [2]
-            assert(self.vq[0]=='pressure')
-            self.switch_V, self.cname, self.vname = 0, 'Q', 'p'
-        elif self.cq[0] == 'pressure':
-            self.c_ids = [0]
-            if self.vq[0] == 'flux':
-                self.switch_V, self.cname, self.vname = 0, 'p', 'Q'
-            elif self.vq[0] == 'volume':
-                self.switch_V, self.cname, self.vname = 1, 'p', 'V'
+        self.v_ids, self.c_ids = [], []
+        self.switch_V, self.cname, self.vname = [], [], []
+
+        for n in range(self.num_models):
+
+            self.v_ids.append(4*n+0)
+            if self.cq[n] == 'volume':
+                self.c_ids.append(4*n+2)
+                assert(self.vq[0]=='pressure')
+                self.switch_V.append(1), self.cname.append('V'), self.vname.append('p')
+            elif self.cq[n] == 'flux':
+                self.c_ids.append(4*n+2)
+                assert(self.vq[n]=='pressure')
+                self.switch_V.append(0), self.cname.append('Q'), self.vname.append('p')
+            elif self.cq[n] == 'pressure':
+                self.c_ids.append(n)
+                if self.vq[n] == 'flux':
+                    self.switch_V.append(0), self.cname.append('p'), self.vname.append('Q')
+                elif self.vq[n] == 'volume':
+                    self.switch_V.append(1), self.cname.append('p'), self.vname.append('V')
+                else:
+                    raise ValueError("Unknown variable quantity!")
             else:
-                raise ValueError("Unknown variable quantity!")
-        else:
-            raise NameError("Unknown coupling quantity!")
+                raise NameError("Unknown coupling quantity!")
 
         self.set_solve_arrays()
 
 
     def equation_map(self):
 
-        self.varmap = {self.vname : 0, 'g' : 1, 'q' : 2, 's' : 3}
-        self.auxmap = {self.cname : 2}
+        self.varmap, self.auxmap = {}, {}
+        for n in range(self.num_models):
+            # self.varmap = {self.vname : 0, 'g' : 1, 'q' : 2, 's' : 3}
+            # self.auxmap = {self.cname : 2}
+            self.varmap[self.vname[n]] = 4*n+0
+            self.varmap['g'] = 4*n+1
+            self.varmap['q'] = 4*n+2
+            self.varmap['s'] = 4*n+3
+            self.auxmap[self.cname[n]] = 4*n+2
 
         self.t_ = sp.Symbol('t_')
-        p_ = sp.Symbol('p_')
-        g_ = sp.Symbol('g_')
-        q_ = sp.Symbol('q_')
-        s_ = sp.Symbol('s_')
-        VQ_ = sp.Symbol('VQ_')
+
+        p_, g_, q_, s_, VQ_ = [], [], [], [], []
+        for n in range(self.num_models):
+            p_.append(sp.Symbol('p_'+str(n+1)))
+            g_.append(sp.Symbol('g_'+str(n+1)))
+            q_.append(sp.Symbol('q_'+str(n+1)))
+            s_.append(sp.Symbol('s_'+str(n+1)))
+            VQ_.append(sp.Symbol('VQ_'+str(n+1)))
 
         # dofs to differentiate w.r.t.
-        self.x_[0] = p_
-        self.x_[1] = g_
-        self.x_[2] = q_
-        self.x_[3] = s_
-        # coupling variables
-        self.c_.append(VQ_)
-        if self.cq[0] == 'pressure': # switch Q <--> p for pressure coupling
-            self.x_[0] = VQ_
-            self.c_[0] = p_
-
-        # df part of rhs contribution (df - df_old)/dt
-        self.df_[0] = self.L*self.C/self.Z * g_ + self.L*self.C * s_
-        self.df_[1] = p_
-        self.df_[2] = VQ_ * self.switch_V
-        self.df_[3] = q_
-
-        # f part of rhs contribution theta * f + (1-theta) * f_old
-        self.f_[0] = (self.L/(self.R*self.Z) + self.C) * g_ + (p_-self.p_ref)/self.R + q_ + (self.L/self.R + self.L/self.Z) * s_
-        self.f_[1] = -g_
-        self.f_[2] = -q_ - (1-self.switch_V) * VQ_
-        self.f_[3] = -s_
-
-        # populate auxiliary variable vector
-        self.a_[0] = self.c_[0]
+        for n in range(self.num_models):
+            self.x_[4*n+0] = p_[n]
+            self.x_[4*n+1] = g_[n]
+            self.x_[4*n+2] = q_[n]
+            self.x_[4*n+3] = s_[n]
+            # coupling variables
+            self.c_.append(VQ_[n])
+        for n in range(self.num_models):
+            if self.cq[n] == 'pressure': # switch Q <--> p for pressure coupling
+                self.x_[4*n+0] = VQ_[n]
+                self.c_[n] = p_[n]
+
+        for n in range(self.num_models):
+
+            # df part of rhs contribution (df - df_old)/dt
+            self.df_[4*n+0] = self.L[n]*self.C[n]/self.Z[n] * g_[n] + self.L[n]*self.C[n] * s_[n]
+            self.df_[4*n+1] = p_[n]
+            self.df_[4*n+2] = VQ_[n] * self.switch_V[n]
+            self.df_[4*n+3] = q_[n]
+
+            # f part of rhs contribution theta * f + (1-theta) * f_old
+            self.f_[4*n+0] = (self.L[n]/(self.R[n]*self.Z[n]) + self.C[n]) * g_[n] + (p_[n]-self.p_ref[n])/self.R[n] + q_[n] + (self.L[n]/self.R[n] + self.L[n]/self.Z[n]) * s_[n]
+            self.f_[4*n+1] = -g_[n]
+            self.f_[4*n+2] = -q_[n] - (1-self.switch_V[n]) * VQ_[n]
+            self.f_[4*n+3] = -s_[n]
+
+            # populate auxiliary variable vector
+            self.a_[4*n+0] = self.c_[n]
 
 
     def initialize(self, var, iniparam):
 
-        var[0] = iniparam[self.vname+'_0']
-        var[1] = iniparam['g_0']
-        var[2] = iniparam['q_0']
-        var[3] = iniparam['s_0']
+        for n in range(self.num_models):
+            try: var[4*n+0] = iniparam[self.vname[n]+str(n+1)+'_0']
+            except: var[4*n+0] = iniparam[self.vname[n]+'_0']
+            try: var[4*n+1] = iniparam['g'+str(n+1)+'_0']
+            except: var[4*n+1] = iniparam['g_0']
+            try: var[4*n+2] = iniparam['q'+str(n+1)+'_0']
+            except: var[4*n+2] = iniparam['q_0']
+            try: var[4*n+3] = iniparam['s'+str(n+1)+'_0']
+            except: var[4*n+3] = iniparam['s_0']
 
 
     def initialize_lm(self, var, iniparam):
 
-        if 'p_0' in iniparam.keys(): var[0] = iniparam['p_0']
+        for n in range(self.num_models):
+            if 'p_0' in iniparam.keys(): var[n] = iniparam['p_0']
+            if 'p'+str(n+1)+'_0' in iniparam.keys(): var[n] = iniparam['p'+str(n+1)+'_0']
 
 
     def print_to_screen(self, var, aux):
 
         if self.ode_parallel: var_arr = allgather_vec(var, self.comm)
         else: var_arr = var.array
 
-        utilities.print_status("Output of 0D model (4elwindkesselLpZ):", self.comm)
+        for n in range(self.num_models):
+            utilities.print_status("Output of 0D model (4elwindkesselLpZ) "+str(n+1)+":", self.comm)
 
-        utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format(self.cname,' = ',aux[0]), self.comm)
+            utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format(self.cname[n],' = ',aux[4*n+0]), self.comm)
 
-        utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format(self.vname,' = ',var_arr[0]), self.comm)
-        utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format('g',' = ',var_arr[1]), self.comm)
-        utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format('q',' = ',var_arr[2]), self.comm)
+            utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format(self.vname[n],' = ',var_arr[4*n+0]), self.comm)
+            utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format('g',' = ',var_arr[4*n+1]), self.comm)
+            utilities.print_status('{:<1s}{:<3s}{:<10.3f}'.format('q',' = ',var_arr[4*n+2]), self.comm)