Ballooning Graphics Memory Space in Full GPU Virtualization Environments

Advances in virtualization technology have enabled multiple virtual machines (VMs) to share resources in a physical machine (PM). With the widespread use of graphics-intensive applications, such as two-dimensional (2D) or 3D rendering, many graphics processing unit (GPU) virtualization solutions have been proposed to provide high-performance GPU services in a virtualized environment. Although elasticity is one of the major benefits in this environment, the allocation of GPU memory is still static in the sense that after the GPU memory is allocated to a VM, it is not possible to change the memory size at runtime. This causes underutilization of GPU memory or performance degradation of a GPU application due to the lack of GPU memory when an application requires a large amount of GPU memory. In this paper, we propose a GPU memory ballooning solution called gBalloon that dynamically adjusts the GPU memory size at runtime according to the GPU memory requirement of each VM and the GPU memory sharing overhead. The gBalloon extends the GPU memory size of a VM by detecting performance degradation due to the lack of GPU memory. The gBalloon also reduces the GPU memory size when the overcommitted or underutilized GPU memory of a VM creates additional overhead for the GPU context switch or the CPU load due to GPU memory sharing among the VMs. We implemented the gBalloon by modifying the gVirt, a full GPU virtualization solution for Intel’s integrated GPUs. Benchmarking results show that the gBalloon dynamically adjusts the GPU memory size at runtime, which improves the performance by up to 8% against the gVirt with 384 MB of high global graphics memory and 32% against the gVirt with 1024 MB of high global graphics memory.

Introduction
Running graphics-intensive applications that include threedimensional (3D) visualization and rendering in a virtualized environment creates a new challenge for highperformance graphics processing unit (GPU) virtualization solutions. GPU virtualization is a technique that allows multiple virtual machines (VMs) to share a physical GPU and run high-performance graphics applications with performance guarantees.
is method is easy to implement, but requires modification every time the version of the API library or GPU driver changes. is method also cannot provide all GPU functions. Direct pass-through [12,13] allocates a GPU exclusively to a single VM and allows it to directly use the GPU with no intervention by the hypervisor. is method provides a high performance that is similar to the native environment, but GPU sharing among VMs is impossible. Currently, AWS [14] and Azure [15] provide GPU services to VMs through the direct pass-through method or through the NVIDIA GRID [16] which is a GPU virtualization solution at the hardware level.
To solve the problems of the approaches mentioned above, GPU virtualization solutions at the hypervisor level, such as gVirt [17], GPUvm [18], and VMCG [19], have been proposed. e gVirt is a full GPU virtualization technology for Intel's integrated GPUs. Unlike dedicated or discrete GPUs in which dedicated graphic cards have independent graphics memory, integrated GPUs share a portion of the system RAM for graphics memory (or GPU memory). e original gVirt divides the graphics memory into chunks and allocates them exclusively to each VM. As a result, a single host could create up to only a maximum of three VMs. e gScale [20,21] solved this scalability problem by dividing the graphics memory into multiple small slots and letting VMs share the slots. Private and physical graphics translation tables (GTTs) are used to translate the virtual address generated by each VM into the physical address for the graphics memory. When a VM is scheduled to run, the private GTT entries of the corresponding VM are copied to the physical GTT. Every time entries are updated in the physical GTT, the modified contents are synchronized with the private GTT that each VM has.
However, the GPU memory size, which is set during the initial creation of the VM, cannot be changed dynamically.
is causes the following problems. First, a VM must be restarted to change the allocated GPU memory size. e user must restart the VM to execute the GPU application that requires GPU memory larger than the current setting. As a result, the VM stops, and the service interruption is expectable. Second, the VM that occupies GPU memory larger than necessary can degrade the performance of other CPU and GPU applications running on other VMs. Whenever a GPU context switch occurs, the CPU copies the private GTT entries to a physical GTT, which also takes up the CPU time of other VMs. If the memory utilization is low, this unnecessary copying overhead may degrade the performance of other VMs. In addition, the GPU context switch time increases as the GPU memory size of each VM gets larger [22]. ird, as we reported in a previous study [22], small GPU memory size affects the performance of GPU workload, especially when VMs run with graphics operations for rendering or high-resolution display environments.
Although several studies [23][24][25][26][27][28][29] dynamically adjusted the memory allocation of existing VMs, in these studies, memory was taken from a specific VM and allocated to another VM when the physical memory was insufficient. As these approaches assume an environment in which each VM has an independent virtual address space, it is difficult to apply those techniques directly to the full GPU virtualization environment in which the same virtual GPU memory space is shared.
In this paper, we propose a dynamic GPU memory ballooning scheme called gBalloon which dynamically increases or decreases the GPU memory size allocated to each VM at runtime. e gBalloon detects performance degradation due to the lack of GPU memory and increases the GPU memory size allocated to each VM. In addition, the GPU memory size of each VM can be reduced when the overcommitted or underutilized GPU memory of a VM creates additional overhead for the GPU context switch or the CPU load due to GPU memory sharing among the VMs. We implement the gBalloon using the 2016Q4 version of gVirt. As the gScale's GPU memory-sharing technique is also implemented, the gBalloon can scale up to 15 Linux VMs. Using various CPU and GPU benchmarks, we also show that the gBalloon dynamically adjusts the GPU memory size at runtime and outperforms the gVirt (modified gVirt in which the gScale's features are added) by up to 8% against the gVirt with 384 MB of high global graphics memory and 32% against the gVirt with 1024 MB of high global graphics memory.
Although current gBalloon is mainly targeted at Intel's integrated GPU, its design principle can be applied to other architectures as well. For example, the proposed idea can be easily applied to other integrated GPUs from AMD and Samsung, where the system memory is used as GPU memory. In addition, we believe that other discrete GPUs with dedicated memory such as NVIDIA can benefit from the gBalloon since they also use graphics translation table for address translation. However, special care has to be taken to reduce the memory copy overhead across the system bus if the gBalloon is implemented over discrete GPUs. As the gVirt is open source and the access to the source codes for the NVIDIA driver and runtime is limited, we decided to use the gVirt as a software platform to verify our proposed idea. e rest of the paper is organized as follows. In Section 2, we outline the structure of gVirt and present the motivations behind the design of gBalloon. In Section 3, we explain the design and implementation issues of gBalloon. In Section 4, we evaluate the performance of gBalloon and compare it with that of gVirt using various configurations. In Section 5, we discuss related works, and in Section 6, we conclude with suggestions for future work.

Background and Motivation
In this section, we provide an overview of the gVirt and discuss the motivations for the proposed approach.

Overview of gVirt.
e gVirt is a high-performance GPU virtualization solution that provides mediated pass-through capability [17]. e gVirt allows VMs to directly access resources that have a large effect on performance and make other privileged operations be intervened through a hypervisor. Due to the restriction on the number of simultaneous VMs in the gVirt, we modified the original gVirt over Xen hypervisor (XenGT) and added the gScale's scalability features [20]. roughout this paper, we consider the modified gVirt as gVirt.
In the gVirt, the mediator located in Dom0 emulates the virtual GPU (vGPU) for each VM. e mediator schedules the vGPU in a round-robin fashion for fair scheduling among the VMs. Considering the high cost of the GPU context switch, each vGPU is switched at approximately 16 ms interval, which is a speed at which people cannot recognize an image change [17]. When the time quantum allocated to each vGPU is fully used, the mediator saves the registers and memory information of the corresponding vGPU and restores the GPU state for the next vGPU. Because the gVirt does not support preemption (although NVIDIA starts to provide preemption capability at the hardware level from Pascal architecture, the feature is not exposed to user's control at the time of this writing) at the time of this writing, it traces the submitted GPU commands and waits for the results until each vGPU can finish its jobs.
For example, if a vGPU executes GPU kernels longer than 16 ms, it runs without preemption. To prevent each vGPU from overusing the time quantum, the gVirt places a limit on the number of GPU kernels a vGPU is allowed to run within the time quantum.
Intel's global graphics memory is divided into two parts as shown in Figure 1: low global graphics memory and high global graphics memory. Only the GPU can access high global graphics memory, but the CPU and the GPU can access low global graphics memory. e CPU can access low global graphics memory through the aperture mapped to the GPU memory.
us, the amount of low global graphics memory that can be accessed depends on the aperture size. e maximum aperture size currently supported by the motherboard is 512 kB, which is mapped up to 512 MB of the low global graphics memory. Figure 1 shows the memory mapping and management structure between global graphics memory and system memory. In the gVirt (modi ed gVirt in which the gScale's features are added), part of the low global graphics memory is shared by all vGPUs, and the high global graphics memory is divided into 64 MB slots that can also be shared among the vGPUs. e virtual address of the global graphics memory is converted into a physical address through the physical GTT. Each vGPU has a private GTT, which is a copy of the physical GTT corresponding to the allocated low global graphics memory and high global graphics memory. e private GTT entries of the vGPU are synchronized every time each vGPU modi es the physical GTT entries. To activate the vGPU in a GPU context switch, if the entry does not exist in the physical GTT, the state is restored by copying the private GTT to the physical GTT. However, if the vGPU is scheduled out, the CPU cannot access the vGPU through aperture. To solve this problem, the gScale allows the CPU to access the vGPU space at all times through the fence memory space pool, which ensures the proper operation of ladder mapping for mapping the guest physical address to the host physical address [20]. e gVirt framework focuses on the acceleration of graphics-intensive applications such as 2D or 3D rendering over a virtualized environment, rather than general-purpose computing on graphics processing units (GPGPU) computing over clouds.

Motivation.
In current gVirt, 64 MB low global graphics memory and 384 MB high global graphics memory are recommended for Linux VM [20] because those memory sizes are enough to support most GPU workloads without performance degradation and crash from the experiments. However, we showed in a previous study [22] that large high global graphics memory can sometimes increase the performance of GPU workloads. We also observed in the study that large high global graphics memory can increase the possibility of overlapping address spaces, which may incur large overhead in a GPU context switch and thus degrade the performance of other VMs. Furthermore, in an environment where the VMs require large global graphics memory, it is highly likely that their GTT entries do not exist in the physical GTT in a GPU context switch because the GPU memory space is shared among the VMs. If a large amount of GPU memory is allocated although it is not fully utilized, unnecessary copies of the GTT entries can occur in a GPU context switch. is increases the time for the GPU context switch, which also decreases the time for each vGPU to occupy the GPU per unit time. As a result, not only the performance of the GPU applications running on all VMs is degraded but also the time for the CPU to copy the GTT entries increases. erefore, the performance of the CPU applications running on the VMs can be degraded as well.
To con rm this, we conducted two experiments to investigate the e ects of copying GTT entries on the performance of GPU application due to excessive occupation of high global graphics memory. For the two experiments, 384 MB and 1024 MB are used for the high global graphics memory size. It is possible to use other con gurations as long as it is bigger than 384 MB and smaller than the size of physical GPU memory. Furthermore, the size should be multiple of slot size. Figure 2 shows the sum of the frames per second (FPS) for each VM by executing Nexuiz 3D benchmarks from Phoronix Test Suite [30] as we increase the number of VMs from 3 to 15. As shown in Figure 2, when the size of the high global graphics memory is small (384 MB), the VMs start to share the GPU memory from the point when the number of VMs reaches around 10. en, the performance of the VMs degrades as we increase the number of VMs to 12 or 15. However, when the size of the high global graphics memory is large (1024 MB), the VMs start to share GPU memory from the relatively small number of VMs. is means that the copying of the GTT entries in the GPU context switch causes a very large performance degradation. When the number of VMs is 6 or 9, the performance at 1024 MB degraded by approximately 3.5 times compared with the performance at 384 MB.
Overall, the performance of the VMs is highly a ected by the size of the GPU memory allocated to each VM, and the memory size must be adjusted at runtime to optimize the performance.

Design of gBalloon
In this section, we describe the design and implementation of the gBalloon that adjusts the GPU memory size of VMs at runtime. As we identi ed in the previous section, the performance of a GPU application can be degraded due to the static allocation of GPU memory. e gBalloon monitors the lack of GPU memory in VMs and then allocates the required amount of GPU memory to the corresponding VM. Furthermore, the gBalloon also checks the performance of CPU and GPU applications and decreases the GPU memory size when the performance of each VM degrades due to the GPU memory sharing among the VMs. e gBalloon is implemented by modifying the gVirt 2016Q4 release [31]. In the following, we present the GPU memory expansion and reduction strategies implemented in the gBalloon in detail.

GPU Memory Expansion Strategy.
When a VM is created, the GPU driver of the VM obtains the available range of GPU memory from the GPU driver in Dom0 and balloons the requested memory area excluding the space that has been allocated to the VM. en, the GPU driver of the VM searches for memory space excluding the ballooned area when allocating a new memory object. If the space for allocating an object is insu cient, the GPU driver creates an empty space by returning the existing memory objects. As a result, the performance of a GPU application that requires graphics-intensive operations, such as rendering operations, is degraded because the same objects are frequently returned. To reduce this overhead, the gBalloon detects the VMs' lack of GPU memory by tracing the number of memory object returns at runtime and reduces the ballooned area of other VMs for the required amount of memory space so that the VM with the lack of memory can use additional GPU memory. Figure 3 shows the process in which the gBalloon allocates additional GPU memory to the guest. When a guest GPU driver must return an existing object due to the lack of GPU memory, the following four steps take place.
Step 1: the guest requests additional GPU memory space from the host.
Step 2: to expand the GPU memory space with the requested size, the host chooses the optimal strategy that can minimize the GPU memory-sharing overhead based on the GPU memory adjustment algorithm which will be explained later.
Step 3: based on the results of the GPU memory adjustment algorithm, the size of the guest's private GTT is increased.
Step 4: nally, the guest receives information about the GPU memory expansion from the host and shrinks the existing ballooned area so that it can be used to allocate objects. For example, as shown in Figure 3, the high global graphics memory area of vGPU1 is expanded to the right, and the shared memory areas of vGPU1 and vGPU3 are increased accordingly.

GPU Memory Reduction Strategy.
As the GPU memory is expanded by the GPU memory expansion request of a VM, the size of the GPU memory shared among the VMs can also be increased. Consequently, the probability that the entry will not exist in the physical GTT during the GPU context switch is increased, resulting in more GTT entry copies. is degrades the performance of the GPU application. In addition, as the number of entries to be copied is also increased, the CPU consumes more time copying the GTT entries, thus degrading the performance of the CPU application.
e gBalloon monitors the CPU cycles consumed for copying GTTentries during the GPU context switch to check the performance degradation of the GPU application. Pro ling is performed at every t cycles. When N is the number of CPUs in the host and C is the CPU cycles consumed for copying GTT entries, the rate of time R copy consumed by the CPU to copy the GTTentries for unit time t can be expressed as follows: Because one CPU processes the copying of the GTT entries, a total of N CPUs consume the cycles for the unit time. us, the number of CPUs must be re ected in R copy .
To check the performance degradation of the CPU application due to the competition among the vCPUs, the gBalloon uses the steal time. e steal time is the time when the vCPU of a VM exists in the ready queue. Assume that w ij is the steal time of the vCPU j of VM i and s ij is the time when the vCPU j of VM i exists in other queues. en, the rate W of the steal times in all VMs can be expressed as follows: A large value of W means that there is severe competition among the VMs. erefore, the state of a physical machine (PM) can be de ned by the values of W and R copy . If both W and R copy are large, then the performance of the CPU application is being degraded by the copying of the GTT entries. In this case, the host must prevent the performance degradation of the CPU and GPU applications by rejecting the GPU memory expansion requests of the VMs and reducing the GPU memory sharing among the VMs. Whereas, if W is large, but R copy is small, although there is severe competition among the VMs, it is not caused by the copying of the GTT entries. In this case, the performance of the CPU application could be degraded due to the increase in R copy .
us, the GTTsize should be reduced if possible. However, if W is small, but R copy is large, there is no competition among the VMs, but many copies of the GTT entries are occurring. In this case, if the VMs perform CPU applications, then the competition among the vCPUs can become more severe due to the copies of the GTT entries, and thus, the performance of the CPU application can degrade. erefore, the host should try to reduce the sharing of GPU memory as much as possible. Finally, if both W and R copy are small, there is no overhead in the current PM, and the host does not need to take any action. Based on the implications described above, the gBalloon identi es the degree of overhead in the host using the values of W and R copy . f(W, R copy ), which represents the degree of overhead incurred due to the copying of the GTTentries, can be expressed as follows: where α and β are the reciprocals of maximum W and R copy . ose parameters normalize the overall value by giving the same weight to W and R copy . e maximum W and R copy are dependent upon a particular hardware platform and are generally determined through experiments.
e gBalloon calculates f(W, R copy ) every 0.5 seconds and classi es the value with two thresholds Threshold low and Threshold high . Based on the classi cation, the gBalloon applies di erent GPU memory policies. e threshold values are experimentally determined from 0 to 2 as the value of f(W, R copy ) ranges from 0 to 2. For the experiments, we use 0.5 (25% of the maximum value) and 1 (50% of the maximum value) for Threshold low and Threshold high , respectively. If f(W, R copy ) is smaller than Threshold low , the gBalloon approves all requests for GPU memory expansion and does not reduce the GPU memory of the VMs. If f(W, R copy ) is larger than Threshold low , the gBalloon starts to decrease the size of the GPU memory allocated to each VM to reduce the overhead of the GPU memory sharing among the VMs. Instead of reducing the GPU memory of all VMs, the gBalloon reduces the GPU memory of the VM that occupies a large amount of GPU memory and has the lowest GPU memory utilization. If f(W, R copy ) is larger than Threshold high , the gBalloon rejects all requests for GPU memory expansion from the VMs and reduces a space corresponding to two 64 MB slots regardless of the GPU memory usage of each VM. e method for reducing GPU memory is to minimize the GPU memory sharing by using Algorithm 1 described in the next subsection.

GPU Memory Adjustment Algorithm.
e GPU memory is adjusted to minimize the GPU memory sharing with the existing VMs. When the gBalloon decides the number of GPU memory slots and the target vGPU to adjust based on the GPU memory expansion and reduction strategies discussed earlier, the spaces at both sides of the GPU memory space allocated to VMs are increased or decreased.
For example, let us assume that vGPU4 initiated a GPU memory expansion request for two slots when there are four vGPUs from vGPU1 to vGPU4 that require two, two, two, and one slots, respectively. Also assume that there are ve slots, and the GPU memory size of the host is 320 MB as shown in Figure 4. In this case, there are three possible methods for expanding the two slots as requested by vGPU4: (1) expanding two slots to the left side, (2) expanding one slot each to the left and right sides, and (3) expanding two slots to the right side. In the rst case, vGPU4 shares slots 1 and 2 with vGPU1, resulting in two shared slots in total. In the second case, vGPU4 shares slot 2 with vGPU1 and slot 4 with vGPU2 and vGPU3, resulting in three shared slots in total. In the third case, vGPU4 shares slot 4 with vGPU2 and vGPU3 and slot 5 with vGPU3, resulting in three shared slots in total. erefore, a strategy of expanding two slots to the left minimizes the number of shared slots. e gBalloon measures the number of shared slots and expands the GPU memory allocated to the VM to minimize the sharing overhead among the VMs.  L, R: the number of slots to expand/reduce to the left side and the right side Lcount ⟵ the number of occupied slots from P K − i to P K − 1 by V j (7) Rcount ⟵ the number of occupied slots from P K + S K to P K + S K + Q − i − 1 by V j (8) count ⟵ count + Lcount + Rcount (9) end for (10) if minimum > count then (11) minimum ⟵ count; Lcount ⟵ the number of occupied slots from P K to P K + i − 1 by V j (8) Rcount ⟵ the number of occupied slots from P K + S K + i − Q to P K + S K − 1 by V j (9) count ⟵ count + Lcount + Rcount (10) end if (11) end for (12) if maximum < count then (13) maximum ⟵ count (14) L ⟵ i (15) end if (16)  e method for reducing GPU memory is the same as that for expanding it. In Figure 4, when the GPU memory of vGPU3 should be reduced by one slot, reducing slot 4 rather than slot 5 can minimize the GPU memory sharing among the VMs. us, the policy for minimizing the GPU memory sharing can maximize the effect of the predictive-copy technique [21] that copies the GTT entries in advance by predicting the next scheduled vGPU before the GPU context switch. e detailed algorithms for memory expansion and reduction are presented below.

Performance Evaluation
In this section, we compare the performance of the gBalloon with that of the gVirt using various workloads. Table 1 shows the experimental environment for the performance evaluation. e global graphics memory size of the host is set to 4 GB, which consists of 256 MB of low global graphics memory and 3840 MB of high global graphics memory. Dom0 does not share the global graphics memory with other domains, but other guest VMs share 64 MB of the low global graphics memory and 3456 MB of the high global graphics memory excluding the Dom0 area. e low global graphics memory size of every guest VM is set to 64 MB as recommended in [20], whereas the high global graphics memory size is set differently depending on the experiments. e experiments use four 3D benchmarks and four 2D benchmarks. To measure the 3D performance, Lightsmark, Openarena, Nexuiz, and Urbanterror of Phoronix Test Suite [30] and Unigine Valley (valley) [32] that requires many rendering operations are used. To measure the 2D performance, Firefox-asteroids (Firefox-ast), Firefox-scrolling (Firefox-scr), gnome-system-monitor (gnome), and Midori of Cairo-perf-trace [33] are used. e performance of the 3D benchmarks is measured by the average number of FPS, and the performance of the 2D benchmarks is measured by the execution time. Furthermore, the NAS Parallel Benchmark (NPB) [34] is used to measure the overhead for the CPU due to the GPU context switch.

Performance Comparison Using a Single GPU Application.
In this subsection, the performance of the gVirt and the gBalloon is compared when valley that requires many rendering operations is executed by multiple VMs. For the gVirt, the gVirt-384 (a gVirt version with the high global graphics memory set to 384 MB) and the gVirt-1024 (a gVirt version with the high global graphics memory set to 1024 MB) are used for the performance comparison. For the 2D and 3D benchmarks that do not demand many rendering operations, additional GPU memory is not allocated because they require a small amount of high global graphics memory. erefore, valley is used to compare the performance of the dynamic GPU memory expansion policy of the gBalloon with that of the gVirt. To observe the performance variations due to the increase in the number of VMs and the change in the degree of GPU memory sharing, experiments were performed in which the number of VMs was increased by three. Figure 5 depicts the performance of the gVirt-384, the gVirt-1024, and the gBalloon normalized to the performance of the gVirt-384 by the sum of the FPS values of all VMs when there are 3, 6, 9, 12, or 15 VMs. When the number of VMs is six or fewer, the gVirt-1024 shows better performance than the gVirt-384 because the overhead from the GPU memory sharing is small. e gBalloon also shows a similar performance as the gVirt-1024 because the gBalloon allocates the required amount of GPU memory to VMs. In particular, as the number of VMs is increased, the performance of the gVirt-1024 and the gBalloon becomes more similar because the effect of overhead on the performance degradation becomes small, and the effect of performance degradation due to the GPU memory sharing becomes large. For this reason, all implementations show similar performance when the number of VMs is nine or higher.

Performance Comparison Using Multiple GPU
Applications. In this subsection, the performance of gVirt-384, gVirt-1024, and gBalloon is compared by running various types of GPU applications on 15 VMs. As shown in Table 2, 15 VMs run randomly selected 2D and 3D benchmarks. Because various benchmarks are mixed, it is possible to compare the degree of the performance degradation of the GPU applications due to the GPU memorysharing overhead that may occur as the requirements for the GPU memory change. Figure 6 shows the performance comparison when the randomly selected 2D and 3D benchmarks shown in Table 2 are executed by 15 VMs. All performance values are normalized to that of the gVirt-384. Strangely, valley in the gVirt-1024 shows a better performance than the gVirt-384 although the overhead is large due to the GPU memory sharing. However, the performance of the other 3D benchmarks is decreased by 50% or higher and the performance of the 2D benchmarks by 25%.
us, the performance of the total VMs drops by 24%, on average, compared with that of the gVirt-384. In the case of the gBalloon, the performance of valley is guaranteed because the GPU memory size of each VM expands as the required amount of GPU memory increases. Moreover, the gBalloon minimizes the overhead due to the GPU memory sharing by dynamically adjusting the GPU memory size according to the GPU memory usage. As a result, the GPU context switch time is decreased, and the performance of valley is increased by up to 28%. e performance of the other benchmarks is similar to that of the gVirt-384. Figure 7 shows a summary of the performance in all benchmarks. e performance of the gBalloon is higher by 8% than that of the gVirt-384 and higher by 32% than that of the gVirt-1024.

Performance Comparison Using CPU and GPU
Applications. In this subsection, the effect of performance degradation in CPU applications caused by the copying of the private GTT entries is analyzed. Among the 15 VMs, 7 VMs run with CPU workloads, whereas the remaining 8 VMs run with GPU workloads. e CPU workload uses cg Scientific Programming of the NPB benchmark, which is a workload to nd the smallest eigenvalue of the matrix using the conjugate gradient method. Figure 8 shows the performance of the CPU and GPU benchmarks normalized to that of the gVirt-384. In the case of the gVirt-1024, the performance of cg is decreased by 19% compared with that of the gVirt-384. is is because the CPU consumes a great deal of time copying the private GTT entries due to the large amount of GPU memory shared among the VMs. Furthermore, the performance of valley is increased slightly by approximately 4% due to this overhead. In contrast, the gBalloon limits the increase in the GPU memory size of the VMs by detecting the overhead of the CPU and the overhead due to the copying of the GTT entries. As a result, the performance of cg is decreased by approximately 1.5%, and the performance of valley is increased by approximately 2% compared with that of gVirt-384.

Overhead and Sensitivity Analysis.
In this subsection, the performance of the gBalloon and the gVirt in a single VM environment is compared. High global graphics memory sizes of 384 MB and 1024 MB are set for the VMs of the gBalloon and the gVirt, respectively. Figure 9 shows the   adaptive behavior in the GPU memory slots of the VMs over time when valley is performed with the dynamic GPU memory expansion policy of the gBalloon. Valley is composed of 18 scenes in total, and the amount of GPU memory required is di erent for each scene. When valley is executed rst, the required amount of GPU memory is increased sharply, and the gBalloon expands three slots in 2 second intervals. en, one slot is expanded for several scenes. e number of slots is increased up to 13, and the size of the high global graphics memory of the VMs is increased to 832 MB. Figure 10 shows the performance of gBalloon for each scene and the overall performance, which is normalized to that of the gVirt. As shown in Figure 10, the overall performance of the gBalloon is lower by approximately 1.9% than that of the gVirt. is is because the slots are increased one by one, resulting in performance degradation due to the temporary lack of GPU memory despite the sharp increase in the amount of GPU memory required in the rst scene. However, this performance degradation is negligible. From the 10th scene when the number of slots becomes 12, the performance degradation disappears due to the frequent expansion requests and the lack of GPU memory. us, the FPS values of the gVirt and the gBalloon are similar.

Related Works
Kato et al. [35], Wang et al. [36], Ji et al. [37], and Becchi et al. [38] proposed technologies for solving the problem of insu cient GPU memory when compute uni ed device architecture (CUDA) is performed in the NVIDIA GPU environment. When the amount of GPU memory is insu cient, the data in the GPU memory are moved to the system memory to secure space in the GPU memory, which is allocated to the applications. However, this copy operation has large overhead when performed at runtime, and the user must use a modi ed API.
Kehne et al. [39] and Kehne et al. [40] proposed a swap policy for reducing the overhead at runtime and improving the resource fairness among GPU applications and the utilization of GPU memory. GPUswap divides the bu er into chunks of a xed size, randomly selects the chunk of an application that occupies the largest amount of GPU memory, and moves the chunk to the system memory when the GPU memory is insu cient. GPUswap randomly selects chunks from applications that occupy the largest amount of memory. However, because the chunk to be removed from the system memory is randomly selected, the performance of the corresponding applications may be degraded if highly reusable data are removed from the GPU memory. To reduce this overhead, GPrioSwap determines the priorities of the chunks by pro ling the GPU memory access counts of the GPU applications and moves the chunk with the lowest priority when the GPU memory is insu cient.
Studies have also been conducted to prevent program crashes when the GPU is shared between containers. Kang et al. [41] proposed a solution that proposes the amount of GPU memory that can be allocated to each container. When a container asks to use more than the limited GPU memory size, ConVGPU rejects the request. In contrast, when the GPU memory is insu cient, ConVGPU lets the container wait until GPU memory becomes available even if the requested amount of memory is less than the limited memory size. However, these studies were targeted at discrete GPUs whose data are transferred through PCIe bus and cannot be directly applied to the heterogeneous system architecture. In this architecture, the system memory is used as the GPU memory, and the data copying between the CPU and the GPU is carried out through a zero-copy bu er.   To improve the memory efficiency in the hypervisor environment, the memory overcommitment technique that decreases or increases the memory allocated to VMs is used. Waldspurger [23], Zhou et al. [24], Zhao et al. [25], Guo [26], Kim et al. [27], and Lu and Shen [28] periodically profiled the access frequencies of pages by nullifying the translation look-aside buffer (TLB) entries of randomly selected pages. Based on this, the least accessed pages are returned when the amount of memory is insufficient. However, this method has a problem because performance degradation may occur due to the nullification of the TLB entries. Furthermore, in [26], [27], [28], and [29], the problem of an inability to respond to sudden changes in VMs' memory demands due to the cyclic overcommitment exists. To solve this problem, the memory pressure aware (MPA) ballooning [27] applies different memory return policies by distinguishing the degree of memory pressure between an anonymous page and a file page. e MPA reduces the performance degradation caused by page return by setting the page with a high probability of becoming the least accessed as the object of return using Linux active and inactive lists. Furthermore, the MPA responds to the VMs' unexpected memory requests by immediately reallocating the memory to sudden memory requests and returning the memory slowly.
Recently, Park et al. [22] proposed a dynamic memory management technique for Intel's integrated GPU called DymGPU that provides two memory allocation policies: size-and utilization-based algorithms. Although DymGPU improves the performance of VMs by minimizing the overlap of the graphics memory space among VMs and thus reduces the context switch overhead, DymGPU's allocation is still static, and the memory size cannot be changed at runtime.

Conclusion and Future Works
In GPU virtualization, due to the static allocation of GPU memory, the performance of VMs that require more GPU memory can be degraded or the GPU application can crash. e gBalloon, proposed in this paper, improves the performance of VMs due to the lack of GPU memory by dynamically adjusting the GPU memory size allocated to each VM. Moreover, the gBalloon detects the increase in overhead due to GPU memory sharing and reduces the GPU memory size of the VMs that unnecessarily occupy a large amount of GPU memory. Consequently, the GPU context switch time is decreased, and the performance of the GPU applications is increased. Furthermore, the performance of the CPU applications is also guaranteed because the CPU load is reduced. e study demonstrated through experiments that the performance of the gBalloon improved by up to 32% when compared with the performance of gVirt with 1024 MB of high global graphics memory.
Currently, the gBalloon increases or decreases only the spaces at both sides to adjust the GPU memory space allocated to VMs. is problem can be solved by allocating non-consecutive spaces of small slot units rather than consecutive GPU memory spaces to VMs. We are currently investigating this issue.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.


ViTa-SLAM: A Bio-inspired Visuo-Tactile SLAM for Navigation while Interacting with Aliased Environments

RatSLAM is a rat hippocampus-inspired visual Simultaneous Localization and Mapping (SLAM) framework capable of generating semi-metric topological representations of indoor and outdoor environments. Whisker-RatSLAM is a 6D extension of the RatSLAM and primarily focuses on object recognition by generating point clouds of objects based on whisking information. This paper introduces a novel extension to both former works that is referred to as ViTa-SLAM that harnesses both vision and tactile information for performing SLAM. This not only allows the robot to perform natural interaction with the environment whilst navigating, as is normally seen in nature, but also provides a mechanism to fuse non-unique tactile and unique visual data. Compared to the former works, our approach can handle ambiguous scenes in which one sensor alone is not capable of identifying false-positive loop-closures.

I. INTRODUCTION
Robots are often equipped with off-the-shelf sensors like cameras which are use for navigation, however, vision is sensitive to extremes in lighting conditions such as shadows or unpredictable changes in intensity as shown in Fig. 1. Whilst other on-board sensors like laser range finders can be used in such situations they too are impaired by reflective and absorbing surfaces. Situations of sensory paucity also occur in nature, e.g., a rat moving through a maze in ill-lit conditions as illustrated in Fig. 2, however, evolutionary processes has equipped animals to gracefully accommodate such scenarios. Given the coarse vision and visually aliased narrow tunnels through which a rat often needs to navigate, they are known to rely on tactile feedback derived from whiskers aside from vision to decipher their own location [1]. Conventional robots lack such a robust capability to interact with their environment through contact. Thus, biomimetic robots are gaining traction [2] which has now made it possible to harness visual and tactile sensory modalities for informed decision making. However, it still remains unclear how to best process and combine information from disparate sensory modalities to aid in spatial navigation.
Previous works on visuo-tactile sensor fusion, such as [3], address the challenge of navigating in unstructured scenes whilst interacting with objects. For this, a low cost camera is mounted such that its field-of-view overlaps with that of a deformable foam rod to acquire coherent data from two sensory modalities.  Alt et al., then demonstrated the application of such systems in household settings [4]. Bhattacharjee et al., tackle the problem of generating a dense tactile map of the environment in which a humanoid robot is operating [5]. For this, a tactile sleeve is utilized along with a complementary Kinect RGB-D camera with redundant field-of-view. In [6], a biologically-inspired visual scene representation inference mechanism is devised that allows a robot to perform haptic interactions with the objects. Like this, the humanoid robot (primary test platform for this work) can poke, push, grasp and visually observe the objects just like an infant. It was argued that this mechanism gets rid of the need to learn visual scene representations by deploying deep-learning methods to train on passive observations and also highlights the need for environmental interactions which is one of the primary themes of this work. Robotic grasping is a challenging task with researchers typically focusing on purely vision, which provides coarse-grained information about object location, shape, and size, or purely on haptic-feedback which provides complementary fine-grained cues about contact forces, textures, local shape around contact points, and deformability, all of which are critical for evaluating an ongoing grasp. In [7], the authors proposed harnessing both visual and tactile data in an end-to-end framework to evaluate the current grasp and maintain estimates of the required grasp adjustments in the future.
Most of the aforementioned works have a requirement for a redundant field-of-view to obtain coherent information from sensory modalities of varying sensing range. Additionally, these methods are mostly concerned with developing a dense haptic object/scene map using a tactile sensor with complementary vision sensor. Whilst these methods allow for environmental interactions they are primarily designed for tactile object exploration and not for performing spatial navigation using mobile robot platforms know as Simultaneous Localization and Mapping (SLAM). To this end, our preliminary results were presented in [8] which provided promising results when jointly considering visuo-tactile sensory modalities for addressing the SLAM problem. Subsequently, a robust sensor fusion algorithm has been developed to integrate information from unique and non-unique sensory modalities such as cameras and whiskers, respectively. Additionally, performance metrics are presented herewith to compare and evaluate model performance against vision or tactile only sensing.

II. BIO-INSPIRED SLAM
This work draws inspiration from two well-known bioinspired SLAM frameworks: RatSLAM, a rat hippocampal model based visual SLAM architecture [9]; and Whisker-RatSLAM, an extension of RatSLAM aimed primarily at tactile object exploration [10]. This work relies on slightly modified variants of these models that are summarized here, followed by a description of the proposed ViTa-SLAM.
For each of the models, an overall system architecture is provided using the following convention: nodes represented by right isosceles triangles represent raw sensory data; nodes represented by ellipse(s) represent pre-processing of sensory data before they are converted to input features represented by rounded boxes. The outputs from the models are represented by regular boxes, the pre-processing and feature generation stages are highlighted in light blue and light red blocks, respectively.

A. Visual-SLAM
When investigating the way rodents navigate from a bioinspired perspective, RatSLAM as introduced in [9], [11]- [14], has been proven to be a capable visual SLAM method. RatSLAM is loosely based on the neural processes underlying navigation in the rodent (primarily rat) brain, more specifically the hippocampus. Fig. 3 shows an overview of the visual-SLAM implementation used in this work.
During the preprocessing phase, the input of a camera (visual data) is downsampled to reduce computational cost and to simulate the coarse vision of rats. In this process the incoming visual data is cropped to remove areas that do not provide unique features, like for example the ground. The cropped image is subsampled and converted to greyscale as shown in Fig. 4.
The preprocessed sensory information is now parsed through 3 major components of the RatSLAM architecture:  The pose cells [9] encode the robot's current best pose estimate. Pose cells are represented by a Continuous Attractor Network (CAN) [11,Ch. 4], the posecell network, to resemble the grid cells as introduced in [15]. The grid cells are neurons found in many mammals and are shown to be used in navigation. In the 3D posecell network, the robot's pose estimate (x, y position and heading angle γ) is encoded as a energy packet that is moved through energy injection based on odometry and place recognition.

Preprocessing Templates
The local view (LV) cells are an expandable array of units used to store the distinct visual scenes as a visual template in the environment using a low resolution subsampled frame of pixels. The visual template generated from the current view is compared to all existing view templates by shifting them relative to each other. If the current view is novel, a local view cell is linked with the centroid of the dominant activity package in the pose cells at the time when a scene is observed. When a scene is seen again, the local view cell injects activity into the pose cells.
The experience map is a semi-metric topological representation of the robot's path in the environment generated by combining information from the pose cells and local view cells into experiences. Each experience is related to the pose cell and local view cell networks via the following 4-tuple: < x, y, γ, V > where x, y, γ represent the location of the cell in the PC network while V corresponds to the view associated with the LV cell that relates to the queried experience [16].
Initially the robot relies on odometry which is subject to an accumulating error. When loop closure events happen, meaning a scene has been seen already, the pose estimate based on the odometry is compared to the pose of the experience and graph relaxation is applied [16].
The following differences to RatSLAM have been introduced in the visual-SLAM implementation. First, we use odometry from the robot instead of visual odometry as was originally done to determine the translational and rotational speeds of the robot. Second, the method of template matching and generation has been modified to account for multiple sensory modalities. Third, the posecell network (PC) is now capable of handling a wider range of robot motion such as moving sideways.

B. Tactile-SLAM
Whisker-RatSLAM is a 6D tactile SLAM algorithm inspired by RatSLAM. Instead of input from a camera, it uses a tactile whisker-array mounted on a robot as its only sensor [17], [18]. The whisker-array consists of 6 × 4 whiskers, each capable of measuring the point of whisker contact in 3D space, and the 2D deflection force at their base [19]. Whisker-RatSLAM [10] has been demonstrated mapping objects and localizing the whisker array relative to the surface of an object. Similar to RatSLAM, Whisker-RatSLAM generates a semi-metric topological map, the object exploration map, which contains complex 6DOF experience nodes.
In [10], the authors proposed combining these object exploration maps with simple 3DOF experience map generated using RatSLAM with the whisker-input resulting in a topological terrain exploration map with two different types of experience nodes. Fig. 5 shows an overview of the Whisker-RatSLAM algorithm. The tactile data acquired by whisking encompasses 3D contact point cloud of the object (3D Cts.) and the deflection data (Defl.). The point cloud is used to generate the Point Feature Histogram (PFH) while the deflection data is used to generate Slope Distribution Array (SDA). Both PFH and SDA are then fused to obtain a 6D Feature Cell (FC). Similar to the RatSLAM experience map, the pose grid cells and FC that were active at a specific 6D pose of the whisker-array are associated with each other and combined into experience nodes. The experience in this case is defined as the 7-tuple: < x, y, z, α, β, γ, F > where x, y, z, α, β, γ represents the 6D pose including euler angles for orientation and F ← {P F H ∪ SDA} represents the features associated with that experience. The experience node form the object exploration map (Obj. Expl. Map). In order to adapt the activation of the pose cell in accord with the robot motion, the odometry information is also used in the pose grid.
The tactile-SLAM implementation is based on Whisker-RatSLAM, but instead of a 6D posecell network this work uses the same 3D posecell network as the visual-SLAM implementation to allow compatibility and to reduce computation cost for navigation in 3D space. Furthermore, the tactile-SLAM implementation used in this work does not use feature cells, but instead combines the SDA and PFH data into 3D tactile template that are used in a similar way as 3D visual templates. Fig. 5 shows an overview of the tactile-SLAM algorithm. The tactile data acquired by whisking encompasses 3D contact point cloud of the object (3D Cts.) and the deflection data (Defl.). The point cloud is used to generate the Point Feature Histogram (PFH) while the deflection data is used to generate Slope Distribution Array (SDA). Both PFH and SDA are then fused to obtain a tactile template (T). Similar to the RatSLAM experience map, the pose grid cells and T that were active at a specific pose of the whisker-array, are associated with each other and combined into experience nodes. The experiences are, opposed to the 7-tuple used in Whisker-RatSLAM, defined as a 4-tuple: < x, y, γ, T > and T ← {P F H ∪SDA} represents the tactile template associated with that experience. The experience nodes also form a semi-metric experience map similar to the visual-SLAM method. Similar to the visual-SLAM method, the robot's odometry information is also used to move the pose grid.
To generate tactile information using whiskers, one challenge is how to control the whisker-array in order to improve the quality of the sensory information. Previous research on rats [20] has identified a number of whisking strategies that rodents use to potentially improve the sensory information they obtain. One of these strategies is called Rapid Cessation of Protraction (RCP) and refers to the rapid reduction in motor drive applied to the whisker when it makes contact with an object during the protraction phase of exploratory whisking [21]. This effectively reduces the magnitude of bend of the whisker upon contact which in artificial arrays, such as shown in [22], improves the quality of sensory information by constraining the range of sensory response to a region best suited for signal processing. Furthermore, damage to the whiskers from contact is significantly reduced.

C. ViTa-SLAM
In this section, we present the details of our novel visuotactile SLAM algorithm which we refer to as ViTa-SLAM.
The overall system architecture for ViTa-SLAM is shown in Fig. 6: 3 kinds of raw sensory data: tactile, visual, and odometry are now utilized simultaneously. Tactile and visual data are converted into visuo-tactile templates (T, V), respectively and hence, need to be pre-processed. A 3D pose cell network is maintained. The experience in this approach is now defined as a 5-tuple: < x, y, γ, V, T > where V is a visual template and T is a tactile template at the 3D pose given by x, y and γ.
The experience map in this case will be referred to as vita map. In contrast with the conventional experience map, the vita map's nodes contains visual and tactile data. The nodes are termed sparse node if the tactile data is empty and dense node otherwise. As an example, when the whiskers do not make contact, the whisker tactile information is not providing any information while the camera can still acquire novel scene information. When the whiskers are whisking a wall/landmark, both the camera and whiskers yield features that allow the creation of informative dense nodes which greatly help visuotactile SLAM.  Algorithm 1 describes ViTa-SLAM in more detail. The visual and tactile processes are running continuously in parallel to the ViTa-SLAM node, as shown by the -symbol and store the current visual and tactile templates V cur and T cur . The visual process follows the same steps as shown in Fig. 4. The resulting visual template is called V cur .
Similarly, the Tactile process is pre-processing the input data consisting of the xy-deflection angles and the whisker contact points in world frame. In line 2, the point feature histogram (P F H) is generated by creating a N -dimensional histogram of the contact points with B bins per dimension. The resulting histogram is then flattened into a N B histogram. For each whisker, the slope of the xy-deflection between the initial contact to the maximum contact during one whisk cycle 1 is computed. The result is the slope distribution array (SDA). The current tactile template containing P F H and SDA is saved as T cur . The two components can be extracted as [t] P F H and [t] SDA .
In the ViTa-SLAM process, if a whisk cycle has been completed, the data stored in the visual and tactile processes, V cur and T cur , is extracted. In line 4, V cur and T cur are matched against all old visual and tactile templates and the id with the closest match m and the corresponding error are returned. Finally, the error is used to determine if a novel template has been detected or a match with an old template has occurred. In either case, energy is injected at the template with match id m. img ←crop image(img) 3: img ←subsample(img) 4: img ← to greyscale(img) if ≤ τ then 6: inject(match id) , m ← arg min( ) 10: return , m After this process, the current visual and tactile templates are appended to the memory and the template match ID is published for experience map generation. The template matching function computes by comparing the current template to all visual and tactile templates in the memory. The visual error V err is computed as the pairwise sum of absolute differences between V cur and all visual templates in the memory. For the tactile data similarly, the PFH and SDA are treated separately and the respective errors (P F H err and SDA err ) are computed. A weighted sum of all obtained error terms yields the error cur between the current visuo-tactile templates and the ones in the memory as: where, In Eq. (1), |·| L1 represents the L1 norm between the corresponding terms. α, β and γ are scaling factors for the respective errors which represent the standard deviations of the raw sensory data. Finally, the returned match ID (m), is the ID of the combined template with the lowest .

III. EXPERIMENTAL SETUP
In this section, we describe the operational environment and the robot platform that were used for empirical validation of the proposed ViTa-SLAM algorithm.

A. Robot Platform
The robot platform used for this research is called WhiskEye (Fig. 7a) the design of which is based on a previous whiskered robot [18]. WhiskEye is composed of a Festo Robotino body, a 3 DoF neck, and a 3D printed head. Mounted on the head are two monocular cameras and an artificial whisker array consisting of 24 macro-vibrissae whiskers arranged into 4 rows of 6 whiskers distributed around a micro-vibrissae array as the nose. Each macro-vibrissae is instrumented to detect 2D deflections of the whisker shaft measured at its base, and is actuated using a small BLDC motor to reproduce the active whisking behaviour observed in small mammals. Motor commands to each whisker, each DoF in the neck, and to the Robotino body of WhiskEye are issued from a small form factor PC mounted on the body running the ROS execution framework. This allows for candidate control architectures to be developed and deployed using the ROS interface to the physical platform or the Gazebo simulation (shown in Fig. 7b) of WhiskEye as used in this study.

B. Operational Environment
As a proof of concept, the algorithm was primarily tested in simulated aliased environments. For this, a 6 × 6 m 2 arena with 4 wall-mounted visual and 3 tactile landmarks was designed as shown in Fig. 8. Within this environment, the robot was made to revolve around the center of the arena, whilst facing outwards, with a radius of 1 m.

IV. PERFORMANCE METRICS
The following metrics were used to evaluate the performance of the proposed ViTa-SLAM against the Visual-SLAM and Tactile-SLAM.

1) Localization Error Metric (LEM)
The localization error metric (LEM) measures the root mean squared error (RMSE) between the true pose and the estimated pose where the error is calculated separately for position and orientation. Thus, where, In Eq. (2),(·) refers to estimated position (orientation) while (·) refers to true position (orientation).

2) Experience Metric (ExM)
The Experience metric (ExM) introduced in [22] provides a performance measure for algorithms like RatSLAM that produce semi-metric topological maps with loop closures. The ExM is comprised of the average rate of re-localization (ARR) and the average rate of correct re-localization (ARCR). The ARR is defined as the ratio of re-localizations to total number of experiences excluding the base set. The base set is a set of initial experiences that has to be selected to define the main loop closure and to provide a reference for future experience to relocalize with The ARCR is defined as the ratio of correct re-localizations to the total number of re-localizations. High values (close to 1) for both factors indicate high certainty in the pose estimate. In order to determine if a re-localization is correct or incorrect, a threshold is used to compare the accuracy of the estimated pose to the ground truth pose. An experience following an incorrect re-localization is labeled invalid until a correct re-localization occurs. An experience following a correct re-localization is labeled valid.

3) Energy Metric (EnM)
In [25], energy metric was proposed as a generic metric for evalution of variety of SLAM algorithms. The SLAM performance was measured in terms of the energy required to transform the SLAM trajectory to the true trajectory. Let N represent the number of relations between experiences in vita map and their corresponding sample points from the set of collected pose data. Then, δ i,j = x i x j represents the transformation from node x i to x j . If T (·) and R(·) represent the translation and rotation operations, then the energy metric (EnM) can be defined as:

V. PERFORMANCE EVALUATION
In this section, we compare the performance of ViTa-SLAM against Visual and Tactile SLAM approaches using the metrics   described above. The experience maps hence obtained are shown in Fig. 9 while the empirical summary of all the metrics is given in Table I. 2 The Energy Metric (EnM) shows that ViTa-SLAM requires less energy to to transform the trajectory to the ground truth by an order of magnitude. This confirms what can be seen from Fig. 9, the experience maps of visual and tactile only SLAM are highly skewed as the result of wrong loop closures.
To further evaluate the quality of loop closure detection we use the Experience Metric (ExM) with the thresholds for position and angular accuracy set to 0.08 m and 4.6 • with the base set defined as the experiences generated during the first full rotation. The ARR and ARCR for ViTa-SLAM show that ViTa-SLAM is able to re-localize more often and correctly in all cases while the other methods always fail. The reason for this failure becomes evident in Fig. 10a and 10b. For the visual-and tactile-SLAM methods, false positive re-localizations already occur in the base set as a result of the aliased environment. As opposed to this, ViTa-SLAM successfully completes one rotation and correctly closes the loop.
The LEM further confirms these findings as the mean pose estimation accuracy of ViTa-SLAM is clearly lower than the state-of-the-art methods.

VI. CONCLUSION AND FUTURE WORKS
This work demonstrated a novel bio-inspired multi-sensory SLAM mechanism for a robot exploring and interacting with 2 Video demonstration available here an environment that presents ambiguous cues. While previous attempts had been made to propose bio-inspired multi-sensory fusion, no prior research allowed for either environmental interactions through contact or fusion of unique and nonunique sensory information. To this end, ViTa-SLAM was presented which utilizes long-range visual and short-range whisker (tactile) sensory information for efficient bio-inspired SLAM. When comparing against earlier approaches that use only vision like RatSLAM or only tactile information like the Whisker-RatSLAM, it was shown that visuo-tactile sensor fusion can handle ambiguities that would otherwise lead to false positive loop-closure detection.
In the future, we plan to investigate bio-inspired sensor fusion approaches based on predictive coding. We also plan to improve the robustness and acuity of the whisking behaviour whilst incorporating spatial attention mechanisms as is seen in rats [26]. Additionally, active spatial exploration strategies will be explored to improve the accuracy of localization and speed of mapping.


Space–time VMS flow analysis of a turbocharger turbine with isogeometric discretization: computations with time-dependent and steady-inflow representations of the intake/exhaust cycle

Many of the computational challenges encountered in turbocharger-turbine flow analysis have been addressed by an integrated set of space–time (ST) computational methods. The core computational method is the ST variational multiscale (ST-VMS) method. The ST framework provides higher-order accuracy in general, and the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the rotor surface. The ST slip interface (ST-SI) method enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST Isogeometric Analysis enables more accurate representation of the turbine geometry and increased accuracy in the flow solution. The ST/NURBS Mesh Update Method enables exact representation of the mesh rotation. A general-purpose NURBS mesh generation method makes it easier to deal with the complex geometries involved. An SI also provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes, and that is enabling the use of nonmatching NURBS meshes in the computations. The computational analysis needs to cover a full intake/exhaust cycle, which is much longer than the turbine rotation cycle because of high rotation speeds, and the long duration required is an additional computational challenge. As one way of addressing that challenge, we propose here to calculate the turbine efficiency for the intake/exhaust cycle by interpolation from the efficiencies associated with a set of steady-inflow computations at different flow rates. The efficiencies obtained from the computations with time-dependent and steady-inflow representations of the intake/exhaust cycle compare well. This demonstrates that predicting the turbine performance from a set of steady-inflow computations is a good way of addressing the challenge associated with the multiple time scales.

Introduction
The challenges faced in computational flow analysis of turbocharger turbines include unsteady flow through a complex geometry, the need for high-resolution flow representation near the rotor surface, high Reynolds numbers, and multiscale flow behavior. The flow unsteadiness comes from the intake/exhaust cycle and the flow in the turbine. An additional challenge is associated with the multiple time scales involved. The computational analysis needs to cover a full intake/exhaust cycle, which is much longer than the turbine rotation cycle because of high rotation speeds. Therefore computations with time-dependent representation of the intake/exhaust cycle require long-durations in the turbine time scale.
The parts of these challenges faced in turbine computations in a more general context have been addressed also by other researchers, with approaches ranging from using a single blade with spatially-periodic boundary conditions (see, e.g., [1][2][3][4][5][6][7][8]) to "sliding interfaces" (see, e.g., [9][10][11][12]). The precursors of the work presented in this article were reported in [13][14][15][16]. The challenge associated with the multiple time scales was addressed in [16] by going through with the long-duration computation. Here we propose to address that challenge in a different way, by calculating the turbine efficiency from computations with steady-inflow representation of the intake/exhaust cycle.

ST-VMS and ST-SUPS
The deforming-spatial-domain/stabilized ST (DSD/SST) method [26][27][28] was introduced for computation of flows with moving boundaries and interfaces (MBI), including fluid-structure interactions (FSI). In MBI computations the DSD/SST functions as a moving-mesh method. Moving the fluid mechanics mesh to follow an interface enables mesh-resolution control near the interface and, consequently, high-resolution boundary-layer representation near fluidsolid interfaces. The stabilization components of the original DSD/SST are the Streamline-Upwind/Petrov-Galerkin (SUPG) [29] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [26] stabilizations, which are used widely. Because of the SUPG and PSPG components, the original DSD/SST is now called "ST-SUPS." The ST-VMS is the VMS version of the DSD/SST. The VMS components of the ST-VMS are from the residual-based VMS (RBVMS) method [30][31][32][33]. The ST-VMS has two more stabilization terms beyond those in the ST-SUPS, and the additional terms give the method better turbulence modeling features. The ST-SUPS and ST-VMS, because of the higher-order accuracy of the ST framework (see [17,18]), are desirable also in computations without MBI.
In the flow analysis presented here, the ST framework provides higher-order accuracy in a general context. The VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the rotor surface.

ST-SI
The ST-SI was introduced in [20], in the context of incompressible-flow equations, to retain the desirable moving-mesh features of the ST-VMS and ST-SUPS when we have spinning solid surfaces, such as a turbine rotor. The mesh covering the spinning surface spins with it, retaining the high-resolution representation of the boundary layers. The starting point in the development of the ST-SI was the version of the ALE-VMS for computations with sliding interfaces [45,46]. Interface terms similar to those in the ALE-VMS version are added to the ST-VMS to account for the compatibility conditions for the velocity and stress at the SI. That accurately connects the two sides of the solution. An ST-SI version where the SI is between fluid and solid domains was also presented in [20]. The SI in this case is a "fluid-solid SI" rather than a standard "fluidfluid SI" and enables weak enforcement of the Dirichlet boundary conditions for the fluid. The ST-SI introduced in [21] for the coupled incompressible-flow and thermaltransport equations retains the high-resolution representation of the thermo-fluid boundary layers near spinning solid surfaces. These ST-SI methods have been applied to aerodynamic analysis of vertical-axis wind turbines [11,12,20], thermo-fluid analysis of disk brakes [21], flow-driven string dynamics in turbomachinery [117][118][119], flow analysis of turbocharger turbines [13][14][15][16], flow around tires with road contact and deformation [113,[120][121][122], lubrication fluid dynamics [123], aerodynamic analysis of ram-air parachutes [124], and flow analysis of heart valves [110,114,115].
In the ST-SI version presented in [20] the SI is between a thin porous structure and the fluid on its two sides. This enables dealing with the porosity in a fashion consistent with how the standard fluid-fluid SIs are dealt with and how the Dirichlet conditions are enforced weakly with fluid-solid SIs. This version also enables handling thin structures that have T-junctions. This method has been applied to incompressibleflow aerodynamic analysis of ram-air parachutes with fabric porosity [124]. The compressible-flow ST-SI methods were introduced in [125], including the version where the SI is between a thin porous structure and the fluid on its two sides. Compressible-flow porosity models were also introduced in [125]. These, together with the compressible-flow ST SUPG method [127], extended the ST computational analysis range to compressible-flow aerodynamics of parachutes with fabric and geometric porosities. That enabled ST computational flow analysis of the Orion spacecraft drogue parachute in the compressible-flow regime [125,126].
In the computations here, the ST-SI enables spinning the mesh covering the solid surface together with the surface. This retains the high-resolution representation of the flow near the surface, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution.

ST-IGA and STNMUM
The ST-IGA is the integration of the ST framework with isogeometric discretization. It was introduced in [17]. Computations with the ST-VMS and ST-IGA were first reported in [17] in a 2D context, with IGA basis functions in space for flow past an airfoil, and in both space and time for the advection equation. Using higher-order basis functions in time enables getting full benefit out of using higher-order basis functions in space. This was demonstrated with the stability and accuracy analysis given in [17] for the advection equation.
The ST-IGA with IGA basis functions in time enables a more accurate representation of the motion of the solid surfaces and a mesh motion consistent with that. This was pointed out in [17,18] and demonstrated in [22][23][24]. It also enables more efficient temporal representation of the motion and deformation of the volume meshes, and more efficient remeshing. These motivated the development of the STN-MUM [22][23][24]; the name was given in [25]. The STNMUM has a wide scope that includes spinning solid surfaces. With the spinning motion represented by quadratic NURBS in time, and with sufficient number of temporal patches for a full rotation, the circular paths are represented exactly. A "secondary mapping" [17,18,22,37] enables also specifying a constant angular velocity for invariant speeds along the circular paths. The ST framework and NURBS in time also enable, with the "ST-C" method, extracting a continuous representation from the computed data and, in large-scale computations, efficient data compression [19,21,113,[117][118][119]128]. The STNMUM and the ST-IGA with IGA basis functions in time have been used in many 3D computations. The classes of problems solved are flapping-wing aerodynamics for an actual locust [22,23,37,101], bioinspired MAVs [24,99,100,102] and wing-clapping [103,104], separation aerodynamics of spacecraft [93], aerodynamics of horizontal-axis [10,25,99,100] and vertical-axis [11,12,20] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [19,113], thermo-fluid analysis of disk brakes [21], flow-driven string dynamics in turbomachinery [117][118][119], and flow analysis of turbocharger turbines [13][14][15][16].
The ST-IGA with IGA basis functions in space enables more accurate representation of the geometry and increased accuracy in the flow solution. It accomplishes that with fewer control points, and consequently with larger effective element sizes. That in turn enables using larger time-step sizes while keeping the Courant number at a desirable level for good accuracy. It has been used in ST computational flow analysis of turbocharger turbines [13][14][15][16], flow-driven string dynamics in turbomachinery [118,119], ram-air parachutes [124], spacecraft parachutes [126], aortas [110,111], heart valves [110,114,115], tires with road contact and deformation [121,122], and lubrication fluid dynamics [123]. The image-based arterial geometries used in patient-specific arterial FSI computations do not come from the zero-stress state (ZSS) of the artery. A number of methods [35,37,[129][130][131][132][133][134][135][136][137][138] have been proposed for estimating the ZSS required in the computations. Using IGA basis functions in space is now a key part of some of the newest ZSS estimation methods [136][137][138][139] and related shell analysis [140].
In the flow analysis presented here, the ST-IGA enables more accurate representation of the turbine geometry, increased accuracy in the flow solution, and using larger timestep sizes. The STNMUM enables exact representation of the mesh rotation.

General-purpose NURBS mesh generation method
To make the ST-IGA use, and in a wider context the IGA use, even more practical in computational flow analysis with complex geometries, NURBS volume mesh generation needs to be easier and more automated. To that end, a generalpurpose NURBS mesh generation method was introduced in [14]. The method is based on multi-block-structured mesh generation with existing techniques, projection of that mesh to a NURBS mesh made of patches that correspond to the blocks, and recovery of the original model surfaces. The recovery of the original surfaces is to the extent they are suitable for accurate and robust fluid mechanics computations. The method is expected to retain the refinement distribution and element quality of the multi-block-structured mesh that we start with. Because there are ample good techniques and software for generating multi-block-structured meshes, the method makes general-purpose mesh generation relatively easy. Mesh-quality performance studies for 2D and 3D meshes, including those for complex models, were presented in [15]. A test computation for a turbocharger turbine and exhaust manifold was also presented in [15], with a more detailed computation in [16]. The mesh generation method was used also in the pump-flow analysis part of the flowdriven string dynamics presented in [119]. The performance studies and test computations demonstrated that the generalpurpose NURBS mesh generation method makes the IGA use in fluid mechanics computations even more practical. The general-purpose NURBS mesh generation method is used also in the turbocharger-turbine flow analysis presented here.

ST-SI-IGA
An SI provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes. This type of mesh generation flexibility is especially valuable in complex-geometry flow computations with isogeometric discretization, removing the matching requirement between the NURBS patches without loss of accuracy. This feature was used in the flow analysis of heart valves [110,114,115], turbocharger turbines [13][14][15][16], spacecraft parachute compressible-flow analysis [126], and flow around tires with road contact and deformation [122]. It is used also in the turbocharger-turbine flow analysis presented here.

Steady-inflow representation of the intake/exhaust cycle
In calculating the turbine efficiency from the computations with the steady-inflow representation of the intake/exhaust cycle, we first perform a set of steady-inflow computations at different flow rates. We then obtain the turbine efficiencies associated with those computations and tabulate them as a function of the tip-speed ratio. Next we calculate the tipspeed ratio for the moving-averaged inflow velocity coming from the time-dependent representation of the intake/exhaust cycle. The turbine efficiency corresponding to that timedependent tip-speed ratio is obtained by interpolation from the tabulated values of the steady-inflow turbine efficiencies.

Computations presented
For comparison purposes, we include from [16] some of the results reported for the computation with the time-dependent representation of the intake/exhaust cycle. The steady-inflow computations are carried out at 12 different flow rates.

Outline of the remaining sections
The governing equations are given in Sect. 2. For completeness, in Sect. 3 we describe the ST-VMS. The computations are presented in Sect. 4, and the concluding remarks in Sect. 5. In the Appendix, we describe the ST-SI and give the stabilization parameters used in the ST-VMS and ST-SI.

Governing equations
Let Ω t ⊂ R n sd be the spatial domain with boundary Γ t at time t ∈ (0, T ), where n sd is the number of space dimensions. The subscript t indicates the time-dependence of the domain. The Navier-Stokes equations of incompressible flows are written on Ω t and ∀t ∈ (0, T ) as where ρ, u and f are the density, velocity and body force. The stress tensor σ σ σ (u, p) = − pI + 2με ε ε(u), where p is the pressure, I is the identity tensor, μ = ρν is the viscosity, ν is the kinematic viscosity, and the strain rate ε ε ε(u) = ∇ ∇ ∇u + (∇ ∇ ∇u) T /2. The essential and natural boundary conditions for Eq. (1) are represented as u = g on (Γ t ) g and n · σ σ σ = h on (Γ t ) h , where n is the unit normal vector and g and h are given functions. A divergence-free velocity field u 0 (x) is specified as the initial condition.

ST-VMS
For completeness, we include, mostly from [17], the ST-VMS: where are the residuals of the momentum equation and incompressibility constraint. The test functions associated with the velocity and pressure are w and q. A superscript "h" indicates that the function is coming from a finite-dimensional space. The symbol Q n represents the ST slice between time levels n and n + 1, (P n ) h is the part of the lateral boundary of that slice associated with the stress boundary condition h, and Ω n is the spatial domain at time level n. The superscript "e" is the ST element counter, and n el is the number of ST elements. The functions are discontinuous in time at each time level, and the superscripts "−" and "+" indicate the values of the functions just below and just above the time level.

Remark 1
The ST-SUPS can be obtained from the ST-VMS by dropping the eighth and ninth integrations.
We also use in the computations the ST-SI, which we include in Appendix A. The stabilization parameters, embedded in the ST-VMS and ST-SI, are included in Appendix B.

Problem setup
The computation with the time-dependent representation of the full intake/exhaust cycle was reported in [16]. Here we carry out the computations with the steady-inflow representation of the intake/exhaust cycle, at 12 different values of the tip-speed ratio: r ω/U IN , where r ω and U IN are the tip speed and inflow velocity. We will compare the efficiencies extracted from the computations with these two representations. Figure 1 shows the turbocharger-turbine model used in the computations with the steady-inflow representation. The model consists of the volute, stator, turbine and the exhaust pipe. The computations with the time-dependent representation had an exhaust gas purifier, instead of an exhaust pipe, and it also had an exhaust manifold. When we compare the efficiencies from the computations with the two representations, we use the same flow region.
The rotor diameter is 30 mm and the rotor speed is 30,000 rpm, which translates to a turbine rotation period of T TR = 2.0×10 −3 s. The gas density and kinematic viscosity are 0.9 kg/m 3  We note that for the operational conditions we have here, incompressible-flow modeling is reasonable. Besides, comparing the laboratory test data from scaled-up models to data from actual-scale models, it was observed that the Mach number made little difference in the turbine efficiency

Efficiency definition
In reporting the turbine efficiency, we first define two quantities, one instantaneous and the other interval-based: where Γ INF and Γ OUTF are the inflow and outflow boundaries used in the efficiency calculation. Similarly, where P(t) is the instantaneous power extracted from the turbine. Then, a general form of the interval-based efficiency is written as We note that we selected the locations of Γ INF and Γ OUTF with the objective of having the smallest possible domain in the efficiency calculation so that we minimize the effect of the manifold and gas purifier.   Figure 2 shows the quadratic NURBS control mesh for the turbine model. The mesh was generated with the generalpurpose NURBS mesh generation method [14,15] described in Sect. 1.4. Table 1 shows the number of control points and elements in different parts of the mesh. Figure 3 shows the four SIs of the mesh. Two of the SIs have an actual slip. The other two are just for mesh generation purpose and connect nonmatching meshes. They are used for independent meshing in the volute region of the computational domain. The STNMUM enables exact representation of the mesh rotation.

Computational conditions
In temporal representation of the mesh rotation, we use quadratic NURBS basis functions. There are 90 time steps per rotation, which is equivalent to a time-step size of 2.22×10 −5 s. The number of nonlinear iterations per time step is 4, with 500, 500, 600 and 800 GMRES iterations

Computations with the steady-inflow representation
The boundaries Γ INF and Γ OUTF are shown in Fig. 4. Figure 5 shows, at the 12 tip-speed ratios, the turbine efficiency averaged over one rotation, based on Eq. (9), where t 3 = t 1 , t 4 = t 2 , and t 2 − t 1 = T TR .

Computation with the time-dependent representation
For completeness, we include from [16] some of the results reported for the computation with the time-dependent representation of the intake/exhaust cycle. The symbol T will denote the duration of the intake/exhaust cycle, and T = 6.0×10 −2 s. Figure 6 shows the velocity magnitude for the full cycle. Figure 7 shows the turbine efficiency η IB (0, T , t, t).
The turbine efficiency corresponding to the time-dependent tip-speed ratio is obtained by interpolation from the steadyinflow turbine efficiencies displayed in Fig. 5. We note from Figs. 8, 9 and 10 that if the averaging period is too short or too long, the efficiencies from the steady-inflow and time-dependent representations do not match well. However, if we choose the right averaging period, which is in this case the turbine rotation period, the efficiencies match very well. We also note that the effect of the manifold

Efficiency (%)
Efficiency from time-dependent representation Efficiency from steady-inflow representation and gas purifier is not significant, provided of course that we report the efficiency based on the same flow region, shown in Fig. 4. Overall, the turbine efficiency can be predicted from a set of computations with the steady-inflow representation of the intake/exhaust cycle, with little difference compared to the efficiency from the time-dependent representation.

Concluding remarks
One of the main computational challenges encountered in turbocharger-turbine flow analysis is the multiple time scales involved. The computational analysis needs to cover a full intake/exhaust cycle, which is much longer than the turbine rotation cycle because of high rotation speeds. Therefore computations with time-dependent representation of the intake/exhaust cycle require long-durations in the turbine time scale. In this article we focused on that challenge. Many of the other computational challenges were addressed earlier by an integrated set of ST computational methods, and we used those methods in the computations here. The core computational method is the ST-VMS. The ST framework provides higher-order accuracy in general, the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the unsteady flow, and the moving-mesh feature of the ST framework enables high-resolution computation near the rotor surface. The ST-SI enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST-IGA enables more accurate representation of the turbine geometry and increased accuracy in the flow solution. The STNMUM enables exact representation of the mesh rotation. A general-purpose NURBS mesh generation method makes it easier to deal with the complex geometries involved. An SI also provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes, and that is enabling the use of nonmatching NURBS meshes in the computations.
The challenge associated with the multiple time scales was addressed earlier by going through with the longduration computation. As another way of addressing that challenge, we proposed here to calculate the turbine efficiency from computations with steady-inflow representation of the intake/exhaust cycle. For that, we first performed a set of steady-inflow computations at different flow rates. We then obtained the turbine efficiencies associated with those computations and tabulated them as a function of the tip-speed ratio. Next we calculated the tip-speed ratio for the moving-averaged inflow velocity coming from the timedependent representation of the intake/exhaust cycle. We obtained the turbine efficiency corresponding to that timedependent tip-speed ratio by interpolation from the tabulated values of the steady-inflow turbine efficiencies. The efficiencies obtained from the computations with time-dependent and steady-inflow representations of the intake/exhaust cycle compare well. This demonstrates that predicting the turbine performance from a set of steady-inflow computations is a good way of addressing the challenge associated with the multiple time scales.
Acknowledgements This work was supported in part by Grant-in-Aid for Challenging Exploratory Research 16K13779 from Japan Society for the Promotion of Science; Grant-in-Aid for Scientific Research (S) 26220002 from the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT); Council for Science, Technology and Innovation (CSTI), Cross-Ministerial Strategic Innovation Promotion Program (SIP), "Innovative Combustion Technology" (Funding agency: JST); and Rice-Waseda research agreement. This work was also supported in part by Grant-in-Aid for Early-Career Scientists 19K20287 (first author) and ARO Grant W911NF-17-1-0046 and Top Global University Project of Waseda University (third author). The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. The HPC resources provided by Nagoya University High Performance Computing Research Project for Joint Computational Science also contributed to obtaining the results reported in the paper.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: ST-SI
In describing the ST-SI (see [20]), labels "Side A" and "Side B" will represent the two sides of the SI. The ST-SI version of the formulation given by Eq. (3) includes added boundary terms corresponding to the SI. The boundary terms for the two sides are first added separately, using the test functions w h A and q h A and w h B and q h B . Then, putting together the terms added to each side, the complete set of terms added becomes where and h is the element length at the SI (see Appendix B.4).
Here, (P n ) SI is the SI in the ST domain, v is the mesh velocity, γ = 1, and C is a nondimensional constant. For explanation of the added SI terms, see [20,120].
On solid surfaces where we prefer weak enforcement of the Dirichlet conditions [20,42] for the fluid, we use the ST-SI version where the SI is between the fluid and solid domains. Using the terminology of Sect. 1.2, that is a fluid-solid SI. This version is obtained (see [20]) by starting with the terms added to Side B and replacing the Side A velocity with the velocity g h coming from the solid domain. Then the SI terms added to Eq. (3) to represent the weakly-enforced Dirichlet conditions become

Appendix B: Stabilization parameters
We include from [141] the element metric tensor in space and in the ST framework. These are used in Appendices B.3 and B.4 in calculation of the stabilization parameters and element lengths.

Appendix B.1: Element metric tensor in space
Components of the Jacobian matrix Q are written as where ξ j is the parametric coordinate in jth direction. We first scale it with a matrix D to take into account the polynomial order or other factors such as the dimensions of the element domain in the parametric space: Remark 2 In the case of finite elements with Lagrange polynomials of order p i in ith direction, with a parametric space of − 1 ≤ ξ i ≤ 1, the scaling matrix can be chosen as WithQ and the direction r to which we want to associate the length, we define the element length (see [141]) as where Remark 3 From this derivation, what we get with D = I has been used in many methods for calculating the stabilization parameters (see, for example, [37]). In those methods, a scaling factor taking the polynomial order into account is applied to the element length, and here we do the scaling in the parametric space, for each of the parametric directions.
Sweeping over all the directions represented by r, we obtain the minimum and maximum element lengths: They are equivalent to and where λ max and λ min are the maximum and minimum eigenvalues of the argument matrix.

Remark 4
In the implementation, we take measures to keep the calculated element length between h MIN and h MAX .

Appendix B.2: Element metric tensor in the ST framework
The ST Jacobian matrix is where θ is the parametric coordinate in time, and the mesh The ST scaling matrix is given as and the scaling becomeŝ The ST metric tensor is defined as

Appendix B.4: ST-SI
The element length used in the ST-SI is given as These were introduced in [