arXiv:2609.22814v1 [cs.DC] 19 Sep 2026
From Idle to Urgent: A Resource-Harvested HPC Workflow for High-Fidelity Seismic Estimation Tsuyoshi Ichimura1 , Kohei Fujita1,2 , Hideaki Ito1 , Wataru Sakurai1 , Muneo Hori3 , and Lalith Maddegedara1 1
Earthquake Research Institute and Department of Civil Engineering, The University of Tokyo, Japan 2 RIKEN Center for Computational Science, Japan 3 Japan Agency for Marine-Earth Science and Technology, Japan
Abstract We propose an Urgent Interactive HPC workflow that dynamically integrates high-fidelity 3D nonlinear analysis with surrogate neural networks (NNs) to enable rapid decision-making during large-scale earthquakes. This approach achieves both ”Resource Harvesting”, which utilizes idle computing capacity during non-emergency periods, and immediate response during crises. By developing two specialized HPC kernels, the proposed method reduces energy-to-solution by 76% and improves throughput by 3.7-fold during normal operations to efficiently construct training datasets, while during emergencies, it couples NN-based inverse analysis with physics-based simulations and dynamic refinement to reduce conventional computational costs by over 97.9%, enabling the generation of highly reliable spatial time-history ground motion distributions within 30 minutes post-earthquake.
1
Introduction
Due to natural forces (e.g., river meandering) and artificial modifications (e.g., cutting and filling), complex shallow 3D ground structures with depths on the order of 100 –101 m are often formed near the surface. During large-scale earthquakes, such ground heterogeneity locally amplifies seismic waves, resulting in significantly severe structural damage at specific sites compared to surrounding areas (e.g., [1, 2]). While numerous methods targeting immediate postearthquake emergency response have been proposed (e.g., [3–5]), further enhancing the reliability of rapid post-disaster decision-making in dense modern urban areas requires moving beyond coarse ground motion estimates. It is crucial to advance toward high-fidelity, high-resolution spatial time-history ground
1
motion simulations that physically account for localized 3D subsurface structural characteristics. For instance, considering the Tokyo metropolitan area in Japan (approximately 105 × 105 m), there are 102 –103 scattered locations of concern, with approximately size of 103 × 103 m each, featuring complex local ground structures that require high-fidelity evaluations. Upon the occurrence of a large-scale earthquake, 101 –102 specific sites experiencing strong shaking must be automatically screened from these numerous candidate locations, and detailed spatial ground motion estimations must be completed within a stringent 30-minute urgent decision-making time frame post-event. Recent advancements in dynamic coupling between observation systems and supercomputers via high-speed networks [6], exemplified by ”Edge-to-HPC” and ”Superfacility” concepts, are making it technically feasible to collect real-time ground motion observations from affected areas immediately after an earthquake and stage them directly into HPC environments. However, these collected data consist of sparse observation points on the ground surface, making it difficult to estimate detailed time-history ground motion distributions across the entire region from these data alone. Estimating ground motion distributions from sparse observations constitutes a challenging inverse problem that generally requires nonlinear optimization—repeatedly executing forward 3D nonlinear ground amplification analyses to search for input bedrock motions consistent with observed waveforms. This implies performing computationally expensive 3D nonlinear finite-element dynamic analyses at least 102−3 times per site for numerous target locations immediately after an earthquake, which is prohibitive due to extreme computational costs within severe time constraints, even when employing state-of-the-art flagship HPC resources. On the other hand, while pure end-to-end surrogate neural networks (NNs) that infer damage distributions directly from observed waveforms offer fast execution, they carry risks of physical inconsistency and reduced explainability due to black-box properties, leaving critical reliability concerns in life-safety emergency decision-making. Meanwhile, research on leveraging HPC to overcome computational resource constraints in emergencies and establishing real-time urgent forecasting workflows have been conducted actively in recent years (e.g., [7, 8]). Building upon these prior studies, designing and constructing a computational workflow tailored to the target hardware and system-operational characteristics is expected to drastically reduce the large analytical costs, improve time-to-insight, and enable rapid decision-making support. To address these challenges, we enable both high energy efficiency during normal operations (”Resource Harvesting”) together with immediate responsiveness during crises (”Urgent Computing”), by proposing a computational workflow based on high-fidelity simulations satisfying physics-based equations and integrating rapid surrogate-NN-based inverse estimation and dynamic refinement for drastic reduction in the massive analytical cost. The remainder of this paper is organized as follows. Section 2 presents the overall architecture of the workflow, which combines normal-operation resource harvesting with emergency response and dynamic refinement. Section 3 details the design and performance evaluation of two specialized core kernels developed 2
to mitigate the computational cost of 3D ground amplification analysis, which represents the main bottleneck of the workflow. Section 4 demonstrates the estimation accuracy of the proposed workflow for realistic earthquake scenarios and validates its practical feasibility within an emergency response timeline.
2
Urgent and Interactive Workflow for Quick Earthquake Response
We present an overview of an urgent and interactive computational workflow designed to estimate high-fidelity spatial time-history ground motion distributions within a realistic time frame for immediate post-event damage estimation and decision support. The characteristics of this workflow lie in coupling surrogateNN-based inverse analysis with accuracy-guaranteed, high-cost physics simulations, by selectively utilizing two specialized core kernels (CPUGPU-Harvest and EBEGPU-Urgent), each optimized at the computer architecture level to align with the distinct resource utilization profiles and power constraints during supercomputer operations under normal and emergency response conditions. In conventional frameworks without surrogate NNs, inversely estimating input motions consistent with observed data required iterating computationally heavy forward 3D nonlinear ground amplification analyses 102−3 times for nonlinear optimization within the tight post-event time frame. Such an approach is virtually impossible within practical execution times for real-time emergency forecasting. In contrast, our workflow couples two specialized computational kernels—optimized through high-performance computing techniques—with inverse analysis based on pretrained surrogate NNs. This dramatically reduces the number of costly 3D analyses conducted during the emergency response phase (down to a minimum of a single run), enabling execution within emergency response time frames. Furthermore, by incorporating an on-the-fly dynamic refinement process on residual errors, the workflow delivers rapid, highly reliable insights for disaster prevention and mitigation decision-making. Here we note that direct inference using pure end-to-end surrogate NNs suffers from inherent drawbacks, including potential physical inconsistencies and limited explainability caused by black-box models. In our workflow, these limitations are overcome by feeding the inference results into a physics-based simulation at least once—thereby yielding high-fidelity spatial time-history distributions that strictly and consistently reflect strong nonlinear behaviors, such as complex subsurface structural effects, energy dissipation, and wave trapping— and by embedding a dynamic refinement process. Consequently, this workflow ensures the high level of reliability essential for life-safety emergency decisionmaking. The details of each phase are described below.
2.1
Pre-Earthquake Phase: Resource Harvesting
To generate training datasets for NNs, a vast number of nonlinear ground amplification analyses must be precomputed across numerous target sites. Since this 3
process requires large energy usage, designing for superior energy-to-solution while accommodating fine-grained, opportunistic scheduling that leverages backfill slots in batch schedulers is important. In our workflow, we employ a highly efficient core kernel based on heterogeneous computing (CPUGPU-Harvest), designed to effectively utilize both the compute and memory resources of CPUs and GPUs. This enables automatic and efficient harvesting of idle resources during normal operations across scales ranging from small fragments on a few nodes to large node allocations, thereby establishing an efficient, automated data generation scheme. Here we note that while training datasets must be regenerated in cases such as when the ground structure information is updated, standard job schedulers can be used to allocate idle resources as each of the training data generation jobs are in sizes of several compute nodes for a short elapsed limit time. Using the datasets pre-generated through this process, surrogate NNs that inversely infer input bedrock motions from surface ground motion observations are constructed for each site prior to an event (along with hyperparameter optimization). To enable rapid, on-the-fly fine-tuning (interactive refinement) during the subsequent emergency phase, we adopt a lightweight NN architecture with a constrained parameter count capable of adapting efficiently within a small number of steps (see Appendix for details on the NN architecture).
2.2
Post-Earthquake Phase: Urgent Response and Interactive Refinement
Immediately following an earthquake, observed ground motion data are received from surface monitoring stations, and the following two-stage process is dynamically executed for sites where strong shaking was recorded. Stage 1: Rapid 1st-Order Estimation. Input bedrock motions are instantaneously inferred using the pretrained NNs. Taking these estimated bedrock motions as input, the multi-node GPU parallel kernel (EBEGPU-Urgent)— specifically designed to prioritize the minimization of time-to-solution during emergencies—is immediately executed. By running the high-fidelity 3D nonlinear ground amplification forward analysis only once, high-fidelity spatial ground motion distributions physically reflecting complex subsurface structures and strong soil nonlinearities are consistently estimated and delivered to decisionmakers within the stringent 30-minute decision-making time frame. Stage 2: Interactive Dynamic Refinement. The observed surface ground motions are compared with the reproduced waveforms obtained from Stage 1. Based on user interface interventions by decision-makers (e.g., requesting priority of analysis for specific targeted areas) or automated thresholding of waveform discrepancies (e.g., relative L1 error exceeding a designated threshold), the results of this 1st-order estimation (pairs of input bedrock motions and surface responses) are immediately fed back into the NN training loop as training data to perform dynamic on-the-fly fine-tuning (an additional training step that incorporates observation-backed data while retaining the search capability derived from random wave training sets; see Appendix for details). 4
Bedrock input motions are then re-estimated using the updated NN, and the ground amplification analysis is re-executed. Through this ”interactive adaptive control loop,” predictive accuracy is progressively refined until it aligns with observational data, continuously providing updated estimates to decision-makers.
3
Design and Optimization of HPC Kernels
As the core of the proposed workflow, we present and evaluate the performance of two specialized core kernels: CPUGPU-Harvest, which minimizes energy-tosolution and maximizes throughput during normal operations, and EBEGPU-Urgent, which minimizes latency and time-to-solution during emergencies. Achieving the detailed time-history ground motion distributions targeted by this workflow requires finite-element simulations using unstructured elements sufficiently small to model complex 3D subsurface structures and ensure numerical convergence across frequency ranges critical to structural response. Furthermore, accurately capturing actual nonlinear soil responses requires complex, high-fidelity physical models (constitutive laws), resulting in large computational costs. Prior research to reduce computational costs of 3D nonlinear ground amplification analyses include urban-scale dynamic simulations featured in SC Gordon Bell Prize sessions [9, 10]. However, these designs heavily prioritized time-to-solution based on simplified physical models, exhibiting limitations when optimizing energyto-solution on modern architectures or adapting to the memory-intensive, highfidelity constitutive models addressed in this study. Therefore, developing a new set of core kernels optimized for modern heterogeneous computing environments is essential to realizing the proposed workflow.
3.1
Mathematical Formulation and Computational Bottlenecks
We describe the 3D ground amplification analysis based on [11], which is used as a basis of [9, 10]. Owing to its capability in modeling complex geometries and naturally satisfying stress-free boundary conditions at the ground surface, the nonlinear wave equation—with material properties that evolve over time according to the physical soil model—is discretized via the finite element method [12]. The Newmark-β method [13] is used for time integration, leading to solving the following equation at each time step to evaluate the response of the nonlinear time-evolution problem: 2 n 4 n M + C + K δun = dt2 dt 4 f n − qn−1 + Cn vn−1 + M an−1 + vn−1 , (1) dt
5
with n q = qn−1 + Kn δun , un = un−1 + δun , 2 δun , vn = −vn−1 + dt n 4 n−1 a = −an−1 − dt v + dt42 δun .
(2)
Here, M, Cn , and Kn denote the mass, damping, and stiffness matrices at the n-th time step, respectively. Vectors δun , un , vn , an , f n , and qn represent the nodal vectors of displacement increment, displacement, velocity, acceleration, outer force, and inner force respectively, while dt denotes the time step. Rayleigh damping is employed, where the element damping matrix Cne is expressed using the element mass matrix Me and element stiffness matrix Kne as Cne = αMe + βKne . Here, coefficients α and β are obtained by solving the following leastsquares problem: "Z 2 # fmax α 1 n + 2πf β df , minimize h − 2 2πf fmin where fmax and fmin are the upper and lower target frequencies, respectively, and hn is the damping ratio that evolves according to the physical soil model. Consequently, Eq. (1) requires solving a system of equations involving matrices Kn and Cn , which change at each time step n due to the physical nonlinear soil model. Second-order tetrahedral elements are used as they are suitable for modeling complex geometry and due to the necessity in evaluating strain for the physical model. Semi-infinite absorbing boundary conditions are applied to the bottom and side boundaries. In contrast to the conventional simplified soil constitutive models used in [9–11], this study adopts the multi-spring model [14]—a high-fidelity physical constitutive law—to evaluate complex soil nonlinear responses more faithfully (see Appendix for details). While this significantly enhances physical fidelity, it requires storing and updating a vast number of history data for each element, leading to a memory-intensive and computationally expensive analysis with stringent requirements on memory bandwidth and capacity. Therefore, designing algorithms tailored to modern architectures, as described in the subsequent sections, becomes important.
3.2
CPUGPU-Harvest: Energy-Aware Kernel for Resource Harvesting
A standard baseline implementation for this core kernel would perform multinode GPU computations where history data for the high-fidelity physical model are stored in GPU memory, the target matrix is stored in 3 × 3 block CRS format, and the system is solved using the Conjugate Gradient method with a 3×3 block Jacobi preconditioner (hereafter referred to as CRSGPU; see Algorithm 1). On the other hand, during multi-case runs aimed at generating datasets for
6
Algorithm 1 Baseline method: CRSGPU. D indicates the stiffness matrix evaluated by the multispring method, while θ indicates the nonlinear spring parameters. CRS-PCG indicates 3 × 3 block Jacobi preconditioned conjugate gradient solver with CRS-based matrix-vector products. A and b indicate left hand side matrix and right hand side vector of Eq. (1), respectively. All computation is done in FP64. 1: for n = 1; n ≤ nt ; n = n + 1 do 2: δun ⇐ CRS-PCG(A, bn )@GPU 3: {Dn , θn } ⇐ Multispring(δun , θn−1 )@GPU 4: A ⇐ UpdateCRS(Dn )@GPU 5: end for
NNs, maximizing energy efficiency and total system throughput across the supercomputer takes precedence over minimizing the execution time of a single case. Therefore, extending [15], we developed CPUGPU-Harvest, a multi-node dense kernel that exhibits superior energy-to-solution, leverages CPU resources by taking advantage of recent improvement in CPU-GPU data transfer bandwidths (e.g., [16, 17]), and enables coupled CPU-GPU execution that fits into idle system backfill slots, executing multi-case analyses with minimal resource overhead. Here, CPUGPU-Harvest is expected to be a broadly applicable approach for conducting many cases of memory-intensive simulations, such as those employing high-fidelity constitutive laws. Our point of departure, the method in [15], was specifically designed for a single GH200 node environment equipped with a CPU with large memory and a high-speed GPU (i.e., configuration with a 72-core Grace CPU with 480 GB of memory and an H100 GPU with 96 GB of memory connected via a 900 GB/s NVLink-C2C interconnect), posing challenges in application to general heterogeneous or large-scale environments. In this study, we generalize this memory/compute decoupling scheme for heterogeneous environments with an arbitrary number of nodes, while integrating the Element-by-Element (EBE) method [18] and multigrid preconditioning to achieve high versatility and throughput via concurrent multi-case execution. Task parallel and task splitting approaches are two primary methods for utilizing both compute and memory resources across CPUs and GPUs. In the task parallel approach, dependencies between computational kernels are analyzed to assign independent kernels across CPUs and GPUs, enabling parallel utilization of CPU and GPU resources (e.g., [19]). In contrast, the task splitting approach divides the domain within each kernel into subdomains handled by the CPU and GPU, computing each kernel in parallel (e.g., [20]). See e.g., [21, 22] for related studies on cooperative CPU-GPU computing. In both approaches, tasks are assigned to the CPU or GPU with the compute and memory in a set; consequently, depending on the system, either compute performance or memory capacity becomes a bottleneck, preventing fully effective resource utilization. To address this, [15] decoupled memory and compute, storing large history data in CPU memory while executing computations on
7
high-performance GPUs, thereby leveraging the strengths of both CPU memory and high-speed GPUs (specifically, multi-spring calculations are executed on the GPU while asynchronously sending/receiving history data stored in CPU memory, enabling large-scale computations that exceed GPU memory capacity). Although the method in [15] targeted only systems where CPU memory was significantly larger than GPU memory, this study generalizes the memory/compute decoupling technique to accommodate architectures with arbitrary CPU-to-GPU memory capacity ratios (Algorithm 2). Here, a portion of the history data is stored in CPU memory and the remaining portion in GPU memory, with an adjustable allocation ratio that enables system-wide memory resource utilization regardless of the system’s CPU-to-GPU memory ratio. History data stored in CPU memory is updated by the aforementioned pipelined transfer and GPU computation, whereas that stored in GPU memory is directly updated by the GPU. Furthermore, for sparse matrix-vector multiplication in the solver, we adopt the EBE method, which generates element matrices on the fly and multiplies them with the right-hand-side vector rather than reading the global matrix from memory. Although this increases floating-point operations compared to CRSbased sparse matrix-vector multiplication, reducing the data volume read from memory yields speedup on modern GPUs. Moreover, it eliminates the need to update the CRS matrix at every time step (Algorithm 1, Line 4), providing additional speedups. In addition, eliminating the need to store the global matrix in CRS format frees up memory capacity, allowing multiple analysis cases (m cases) with varying input ground motions to be computed concurrently; this reduces random data accesses during sparse matrix-vector multiplication, leading to further performance gains. Additionally, using multigrid preconditioning in the solver enables efficient solution convergence. By incorporating halo exchanges into sparse matrix-vector multiplication and MPI Allreduce into inner products within the solver, parallel execution across multiple computer nodes is achieved. The developed kernel, CPUGPU-Harvest, keeps high-speed GPUs continuously active while fully utilizing both CPU and GPU memory capacities. Executing such dense computations on a small number of nodes minimizes data access and communication overheads, enabling the efficient harvesting of idle supercomputer capabilities with high energy efficiency.
3.3
EBEGPU-Urgent: Time-to-Solution-Optimized GPU Kernel
During the post-disaster emergency phase, minimizing the latency and time-tosolution for a single analysis case becomes top priority, where slight degradation in energy efficiency can be tolerated. Therefore, by extending CRSGPU, we developed EBEGPU-Urgent, a pure GPU kernel designed to enable urgent decisionmaking within 30 minutes. The specific algorithm of EBEGPU-Urgent (see Algorithm 3) targets a single analysis case and improves throughput by adopting the EBE formulation for sparse matrix-vector multiplication, which also eliminates 8
n } indicates parts of the spring paramAlgorithm 2 CPUGPU-Harvest. {θpart(i) eters allocated on GPU memory (i = 0) or CPU memory (i = 1, 2, ..., np ). Lines 14 and 15 (or lines 17 and 18) are conducted asynchronously such that the CPUGPU transfers and GPU computation are overlapped. ”que” stores the index of multispring data on CPU (excluding those residing on GPUbuffer1 and 2. For example, in the initial step of n = 1, que = {3, 4, ..., np }). EBE-MultiIPCG is a conjugate gradient solver with multi-grid and mixed-precision preconditioner, with element-by-element method used for conducting matrix-vector products. m cases are solved together to reduce random accesses. Computation in preconditioner of EBE-MultiIPCG is done in FP32, while all other computation is done in FP64.
1: // Copy part of CPU data to GPU buffers 0 2: CPU[θpart(1) ]→ GPUbuffer1 0 ]→ GPUbuffer2 3: CPU[θpart(2) 4: for n = 1; n ≤ nt ; n = n + 1 do 5: // Solver 6: δun ⇐ EBE-MultiIPCG(Dn−1 , bn )@GPU 7: // Multispring on GPU memory n−1 n n 8: {Dn part(0) , θpart(0) } ⇐ Multispring(δu , θpart(0) )@GPU
9: // Startup of piplined multispring on CPU memory 10: Compute MS on GPUbuffer1; GPUbuffer1.state ⇐ done 11: // Piplined multispring on CPU memory 12: while que.size() != 0 do 13: if GPUbuffer1.state == done then 0 14: GPUbuffer1→CPU; CPU[θpart(que.pop()) ]→GPUbuffer1 15: Compute MS on GPUbuffer2; GPUbuffer2.state ⇐ done 16: else 0 17: GPUbuffer2→CPU; CPU[θpart(que.pop()) ]→GPUbuffer2 18: Compute MS on GPUbuffer1; GPUbuffer1.state ⇐ done 19: end if 20: end while 21: // Closing of piplined multispring on CPU memory 22: if GPUbuffer1.state == done then 23: Compute MS on GPUbuffer2; GPUbuffer2.state ⇐ done 24: else 25: Compute MS on GPUbuffer1; GPUbuffer1.state ⇐ done 26: end if 27: end for
9
Algorithm 3 EBEGPU-Urgent. D indicates the stiffness matrix evaluated by the multispring method, while θ indicates the nonlinear spring parameters. EBEPCG indicates a 3 × 3 block Jacobi preconditioned conjugate gradient solver with EBE-based matrix-vector products. All computation is done in FP64. 1: for n = 1; n ≤ nt ; n = n + 1 do 2: δun ⇐ EBE-PCG(Dn−1 , bn )@GPU 3: {Dn , θn } ⇐ Multispring(δun , θn−1 )@GPU 4: end for
the need for CRS matrix updates. Furthermore, under the premise of deploying large node allocations, all history data are stored entirely in high-speed GPU memory, and all computations are executed on high-speed GPUs. Overlapping sparse matrix-vector multiplication computation with halo exchange communication further enhances scalability across large node counts. Although CPU cores and memory are not utilized, employing a large number of high-speed GPUs achieves extreme computational speed.
3.4
Performance Evaluation on the Miyabi Supercomputer
The performance of CPUGPU-Harvest and EBEGPU-Urgent is evaluated in comparison with the baseline CRSGPU on Miyabi [23]. Miyabi is a jointly operated supercomputer system of the Joint Center for Advanced High Performance Computing (JCAHPC), managed by the Information Technology Center at The University of Tokyo and the Center for Computational Sciences at the University of Tsukuba. In this study, we utilize Miyabi-G, where each compute node is equipped with a single NVIDIA GH200 Grace Hopper Superchip. Each node has one H100 GPU (96 GB, 4.0 TB/s) and one Grace CPU (120 GB, 512 GB/s) interconnected via NVLink-C2C (900 GB/s). A total of 1,120 compute nodes are connected via InfiniBand NDR200 in a full-bisection Fat-Tree topology. The implementations were performed at the same level of optimization as in [24], which conducts coupled CPU-GPU execution for similar analyses. Measurements reported below were obtained across all time steps (nt = 14500), with energy consumption evaluated using module-wide power values—including the CPU, GPU, and memory subsystems—obtained via nvidia-smi -q -d POWER. First, we compare the performance during resource harvesting in the preearthquake phase. Here, throughput and energy-to-solution are compared when using 4 nodes, which is required by the conventional CRSGPU method to solve the problem within a practical execution time (Table 1). We first compare the baseline CRSGPU with CPUGPU-Harvest executing a single case (m = 1). In CPUGPU-Harvest, half of the multi-spring data is stored in CPU memory and the other half in GPU memory, executing with pipelined overlap between CPU data transfers and GPU computation. Although pipeline execution overhead increases the execution time of the multi-spring (MS) section from 799 s to 879 s, the solver section is reduced from 5,352 s to 3,296 s owing to multigrid preconditioning and the EBE method. Combined with the elimination of CRS
10
Table 1: Performance for Pre-Earthquake Phase measured on 4 nodes. CGH stands for CPUGPU-Harvest. m is the number of cases solved simultaneously. MS indicates time for multispring computation. Time is shown per case.
CRSGPU CGH [m = 1] CGH [m = 4]
Elapsed time in sec. Total (Solver, MS, CRS) 8825 (5352, 799, 2590) 4255 (3296, 879, – ) 2370 (1890, 419, – )
Power /node 630 W 504 W 569 W
Required energy 22.2 MJ 8.58 MJ 5.40 MJ
Table 2: Memory usage for Pre-Earthquake Phase measured on 4 nodes. Per node values obtained by numactl -H are shown. CRSGPU CPUGPU-Harvest [m = 1] CPUGPU-Harvest [m = 4]
CPU memory 19.3 GB 28.4 GB 75.0 GB
GPU memory 43.1 GB 23.4 GB 89.0 GB
matrix updates, the overall application throughput improves by a factor of 2.1 (8825/4255). At this stage, as shown in Table 2, CPUGPU-Harvest (m = 1) leaves room in both CPU and GPU memory, enabling concurrent execution of m = 4 cases (the multi-spring data size is 20.4 GB per node per case, half of which is assigned to the CPU side and half to the GPU side even when m = 4). This reduces random memory accesses, further improving the throughput (execution time per analysis case) of the solver section from 3,296 s to 1,890 s compared to m = 1. For the entire application, this corresponds to an 1.8-fold throughput improvement relative to m = 1, and a 3.7-fold improvement relative to the conventional CRSGPU method (8,825 s → 2,370 s/case). Furthermore, using the EBE method reduces total memory transfers per case relative to CRSGPU, which lowers power consumption from 630 W to 569 W. Consequently, the energy-tosolution is reduced to 0.24 times that of the baseline (a 76% reduction in energy consumption). Note that executing EBEGPU-Urgent on 4 nodes (described below) requires 4,143 s (10.2 MJ), indicating that CPUGPU-Harvest with dense computations delivers superior performance in terms of both throughput and energy-to-solution. Thus, a high-throughput, energy-efficient algorithm tailored to the pre-earthquake phase has been successfully realized. Using this approach, the 100-case training dataset used in the application example in Section 4 can be computed in sets of 4 nodes for 2.6 hours (2, 370 s×4 cases), achieving a granularity that can be flexibly scheduled into short idle backfill slots during normal supercomputer operations. Moreover, the entire dataset can be constructed in a practical resource size of 263 node-hours in total. Next, to evaluate urgent computing performance during the post-earthquake phase, we measure the strong scaling performance of EBEGPU-Urgent (Table 3). Although overhead from strong scaling causes a decrease in parallel efficiency,
11
Table 3: Performance for Post-Earthquake Phase. MS indicates time for multispring computation.
EBEGPUUrgent
Compute nodes 4 8 16 32 64 128
Elapsed time Total (Solver, MS) 4143 s (3295 s, 769 s) 2321 s (1886 s, 389 s) 1411 s (1179 s, 202 s) 965 s (834 s, 111 s) 855 s (774 s, 63 s) 771 s (725 s, 31 s)
Required energy 10.2 MJ 10.7 MJ 10.7 MJ 12.6 MJ 17.5 MJ 24.1 MJ
adhering to the emergency-phase design philosophy of prioritizing latency minimization enables the complete 14,500-step analysis to be computed in just 771 s (12.1 min) when using 128 nodes. The required energy on 128 nodes is 24.1 MJ; thus, EBEGPU-Urgent achieves an 11.4-fold speedup (8825/771) over CRSGPU on 4 nodes (22.2 MJ) while consuming roughly comparable energy, making it wellsuited for immediate post-earthquake deployments where power consumption should be constrained. In this manner, an ultra-fast, single-case time-to-solution tailored to the post-earthquake phase is achieved with high energy efficiency.
4
Validation and Application: High-Fidelity Urban Ground Motion Estimation
In this section, we apply the proposed workflow to immediately estimate spatial ground motion distributions across the target region based on ground motion recorded at a single surface observation point. During large earthquakes, ground motions are typically observed at only a sparse set of surface stations relative to the domain size and spatial resolution required for post-disaster damage estimation; thus, there is strong demand for methodologies capable of rapidly estimating spatial ground motion distributions as demonstrated in this study.
4.1
Simulation Setup and High-Fidelity Geotechnical Model
As an example target site, we consider a location near Yokohama City, Kanagawa Prefecture, Japan. This site features soft sedimentary layers forming a complex 3D structure, resulting in significantly larger ground shaking compared to surrounding areas. That is, body waves are converted into surface waves and trapped, inducing large 3D nonlinear ground amplification, making it an ideal site for the evaluation conducted in this section considering damageinducing 3D nonlinear ground amplification. To understand the ground motion mechanisms at this site and obtain insights for disaster mitigation, a detailed geotechnical structure model was constructed by a dedicated committee through comprehensive investigations using various datasets; thus, this study utilizes this
12
Table 4: Material properties of the soil structure First layer Second layer Bedrock
a)
Vp m/s 700 1400 2100
Vs m/s 100 300 700
ρ kg/m3 1500 1800 2100
Elevation of bedrock (m)
First layer
Second layer
hmax 0.23 0.23 0.01
γr 0.007 0.001 ∞
b)
P2
Bedrock P1
Figure 1: (a) 3D ground structure model, and the position of observation points P1 and P2 . (b) Close-up view of the region shown by the rectangle in (a). ground model. The constructed model forms a three-layer stratified structure over an area of 1,696 m east-west by 1,920 m north-south, with soil material properties summarized in Table 4. In [11], 3D nonlinear analysis results were compared against observation data recorded during the devastating 2011 Tohoku earthquake, demonstrating good validation performance of this model. For this model, similar to [11], frequencies up to 2.5 Hz are targeted to capture the primary damage-inducing components, and a finite-element model using unstructured second-order tetrahedral elements was generated to ensure at least 10 elements (21 nodes) per wavelength (given a minimum shear wave velocity of 100 m/s and a wavelength of 40 m at 2.5 Hz, the minimum element size is set to approximately 4 m based on the 10-element threshold). Figure 1 shows the generated 3D ground finite-element model with complex geometry comprising three soil layers (32,502,492 degrees of freedom and 7,781,075 tetrahedral elements). Here, a Cartesian coordinate system is used with the origin placed at the southwest bottom corner of the model, where the x-, y-, and z-axes represent the east-west, north-south, and vertical directions, respectively. In the 3D nonlinear ground amplification analysis conducted in this study, the time step is set to dt = 0.005 s as in [11], with a relative error convergence threshold of 10−8 , running for 14,500 time steps. To evaluate the estimation performance of the proposed workflow under
13
Table 5: Quantitative accuracy of estimated of input wave Pnt Pnt wave (x-component |yi |, where i is the |xi − yi |/ i=1 and wave at P1 ). Here, Err(x, y) = i=1 time step number and nt is the number of time steps. (a) Kobe scenario
Phase II:
K,(1) Err(ssim,P1 , sK ref,P1 ): 0.128 K,(2) Err(ssim,P1 , sK ref,P1 ): 0.073
Phase I:
Err(ssim,P1 , sC ref,P1 ): 0.119
Phase II:
Err(ssim,P1 , sC ref,P1 ): 0.084
Phase I:
K,(1)
Err(best , bK ref ): 0.110 K,(2)
Err(best , bK ref ): 0.066
(b) Chuetsu-oki scenario C,(1)
Err(best , bC ref ): 0.113
C,(1)
C,(2)
Err(best , bC ref ): 0.081
C,(2)
nonlinear behavior driven by realistic large earthquake motions, and considering the localized domain, uniform plane waves are input from the bottom of the ground model—a practice widely adopted in conventional studies. Specifically, the Kobe wave (hereafter referred to as true input bedrock motion bK ref ) and the Chuetsu-oki wave (hereafter referred to as true input bedrock motion bC ref ) are applied as inputs, considering cases where ground motions are observed at points P1 and P2 shown in Fig. 1. That is, we rapidly estimate spatial ground motion distributions using solely the surface waveform observed at P2 , and validate the estimation accuracy against the waveform observed at P1 . Note that bK ref was generated based on seismic motion recorded at Nakayamate, Chuoku, Kobe City, Hyogo Prefecture during the devastating 1995 Kobe earthquake (Hyogoken-Nanbu earthquake) [25], which is widely used in seismic design in Japan. Meanwhile, bC ref was generated based on seismic motion recorded at Komeda, Izumozaki Town, Niigata Prefecture during the severe 2007 Niigataken Chuetsu-oki earthquake [25]. Specifically, because these records are surface observation data, their amplitudes were halved to convert them into bedrock input motions. Furthermore, following [11], frequency components up to 2.5 Hz were extracted to focus on the frequency band associated with heavy damage (i.e., a 0.2–0.5–2.4–2.5 Hz bandpass filter is applied). bC ref is characterized by having slightly higher frequency content compared to bK ref . Since the waveform evaluation results were similar across x-, y-, and z-components, the x-component is presented as a representative case. Additionally, the refinement process is configured to run whenever the relative L1 error (defined as Err in Table 5) is equal to or exceeds 0.10.
4.2
Case I: Kobe Scenario
First, bK ref was input into the ground model, and a 3D ground amplification analysis using the high-fidelity physical model described in the previous section was conducted (hereafter considered as the reference solution). Figure 2(e) shows the time-history velocity norm at the surface sK ref . We can see that complex 3D nonlinear responses reflecting the complex ground structure occur, exhibiting significant time-history variations and vastly different responses depending on 14
Figure 2: Results for Kobe and Chuetsu-oki scenarios. Kobe scenario: (a) waves at P1 , (b) input wave, and (e) time history response distribution. Chuetsu-oki scenario: (c) waves at P1 , (d) input wave. The x-component is visualized for the waveforms. K the location. From this analysis result, surface waveforms (sK ref,P1 and sref,P2 ) were observed at points P1 and P2 shown in Fig. 1 (Fig. 2(a-1) shows the waveform at P1 ). To examine the effects of the high-fidelity physical model, waveforms computed without enabling the physical model were also observed (Fig. 2(a-2) shows the waveform at P1 ). Specifically, because this physical model induces nonlinearity—resulting in stiffness degradation and damping increase—based on the reference strain γr , setting a large value for γr prevents the physical model from exhibiting nonlinear effects. Therefore, γr = 10 was set here to suppress nonlinearity. In the linear computation result, waves are trapped by the ground structure without undergoing stiffness degradation or damping, leading to significantly severe oscillations. Conversely, under the nonlinear setting, substantial stiffness reduction and damping lead to markedly different waveforms, highlighting the importance of employing a high-fidelity physical model. Phase I: Rapid First-Order Estimation: We attempt immediate ground motion distribution estimation using the proposed workflow. First, we attempt to inversely estimate the input bedrock wave from sK ref,P2 . Although this process is a nonlinear optimization problem to which various methods can be applied, we employ rapid evaluation via a NN surrogate model to achieve immediate post-disaster damage estimation. That is, a NN that estimates input motions
15
at the bottom of the model from surface observation data is constructed and used to infer the input wave. Specifically, 100 random waves were generated with frequency components above 2.5 Hz filtered out and amplitudes uniformly distributed between -0.6 and 0.6 for the x- and y-components and between -0.3 and 0.3 for the z-component (this configuration reflects the fact that vertical components are typically smaller than horizontal components in actual earthquake ground motions). Using these waves as inputs, 3D ground amplification analyses were performed using CPUGPU-Harvest, and response waveforms were recorded at P2 . Through this procedure, 100 pairs of time-history data corresponding to input random waves and their surface responses at P2 were obtained. Using these data pairs, a NN surrogate model was constructed prior to the event using the method described in the Appendix; this model is denoted as N N P2 . Taking K,(1) P2 sK is denoted as best ref,P2 as input, the initial bedrock wave inferred by N N (see Fig. 2(b-2)). Although point P2 is subject to complex subsurface structures K,(1) and strong nonlinearities, N N P2 —which accounts for these effects—infers best K with almost the same envelope and phase properties as that of bref , leading to low error of Err = 0.110 (see Figs. 2(b-1), 2(b-2) and Table 5). Next, a 3D K,(1) ground amplification analysis was executed using the estimated input best . K,(1) Comparing the computed waveform ssim,P1 at the unlearned validation point P1 against the reference solution sK ref,P1 yields Err = 0.128, also demonstrating close agreement (see Figs. 2(a-1), 2(a-3) and Table 5). Phase II: Interactive Refined Update: Next, we attempt to improve the accuracy of the estimated spatial ground motion distribution. Specifically, K,(1) K,(1) newly acquired data (the simulated surface waveform ssim,P2 at P2 and best ) are added to fine-tune N N P2 using the method detailed in the Appendix— adapting the model to observational facts while preserving the search capability P2 derived from the random wave dataset—thereby constructing N NK,(1) . Using K,(2)
P2 with improved accuracy N NK,(1) and sK ref,P2 , an updated bedrock motion best is inversely estimated and used as input to execute a second 3D ground amK,(2) plification analysis. The bedrock motion best and the resulting surface waveK,(2) form ssim,P1 obtained at P1 from this second analysis are shown in Figs. 2(b-3) and 2(a-4). Through this refinement process, the estimation error of the input bedrock motion decreases from 0.110 to 0.066, and the surface response error at the unlearned point P1 also decreases from 0.128 to 0.073, clearly demonstrating an improvement in accuracy (see Table 5). Finally, Fig. 3(a) shows the surface Spectrum Intensity (SI) [26] distributions K,(1) K,(2) obtained when applying bK as inputs (the SI value is a comref , best , and best mon measure used to estimate seismic damage to structures). We can see that the 1st-order estimation gives SI distribution in high accuracy of 0.081/2.7 = 3%, and that the refined update improves the SI prediction in regions with large SI values, which are of particular interest in disaster mitigation. These comparisons indicate that the SI distribution—which directly correlates with damage assessment—can be estimated accurately, and that the refinement works effectively.
16
0.61
SI (m/s) K (a-1) 𝑠𝑠ref
0.68
SI (m/s) C (b-1) 𝑠𝑠ref
0.0 Difference (m/s) 0.081
2.7
K, 1
K (a-2) 𝑆𝑆𝐼𝐼diff 𝑠𝑠sim , 𝑠𝑠ref
(a) Kobe scenario
2.8
0.0
K, 2
K (a-3) 𝑆𝑆𝐼𝐼diff 𝑠𝑠sim , 𝑠𝑠ref
Difference (m/s)
C, 1 C (b-2) 𝑆𝑆𝐼𝐼diff 𝑠𝑠sim , 𝑠𝑠ref
(b) Chuetsu-oki scenario
0.2 C, 2
C (b-3) 𝑆𝑆𝐼𝐼diff 𝑠𝑠sim , 𝑠𝑠ref
Figure 3: SI values obtained p at surface. Here, differences of SI values are computed as SIdiff (sim, ref) = {SIx (sim) − SIx (ref)}2 + {SIy (sim) − SIy (ref)}2 .
17
4.3
Case II: Chuetsu-oki Scenario
Similar to Case I, bC ref was first input into the ground model to perform a 3D nonlinear ground amplification analysis, and the surface reference waveforms C sC ref,P1 and sref,P2 at P1 and P2 were recorded. Because the Chuetsu-oki wave C and Kobe wave possess different characteristics, bK ref and bref are fundamentally distinct (see Figs. 2(b-1) and 2(d-1)). Moreover, even when observed at the C same point P1 , waveforms such as sK ref,P1 and sref,P1 exhibit markedly different properties, highlighting the importance of observing ground motions and conducting rapid post-event damage estimation that accounts for 3D nonlinear responses on an event-by-event basis (see Figs. 2(a-1) and 2(c-1)). Phase I: Rapid First-Order Estimation: We attempt immediate ground motion distribution estimation using the proposed workflow. First, we attempt C to inversely estimate the input bedrock wave from sC ref,P2 . Inputting sref,P2 into the pre-constructed shared surrogate model N N P2 described in Case I, we inC,(1) ferred the bedrock input wave best . As shown in Figs. 2(d-1) and 2(d-2), despite having distinct characteristics from the Kobe wave—specifically having larger high-frequency components—it is accurately estimated by the preconstructed surrogate model trained on random waves (Err = 0.113). Using C,(1) the estimated input best , a 3D ground amplification analysis was executed. C,(1) The computed waveform ssim,P1 at the unlearned observation point P1 also shows good agreement (Figs. 2(c-1) and 2(c-2), Err = 0.119 in Table 5). Phase II: Interactive Refined Update: Next, we attempt to improve the accuracy of the estimated spatial ground motion distribution from the firstC,(1) order estimation. Specifically, using the pair consisting of best and the reproC,(1) duced surface wave ssim,P2 from Phase I, N N P2 is fine-tuned on the fly using the method detailed in the Appendix—adapting the model to observation data while preserving search capability derived from random wave data—to rapidly P2 construct an updated N NC,(1) . Using this updated model, the input bedrock C,(2)
wave best is inversely estimated from the observed wave sC ref,P2 , and used as input to perform the second 3D nonlinear ground amplification analysis. The resulting simulated waveform obtained at P1 from this 3D analysis is denoted C,(2) as ssim,P1 . Comparing these waveforms reveals that discrepancies in the firstorder estimation are refined for both the input wave and the observed waveform at P1 , resulting in improved accuracy (see Figs. 2(d-3) and 2(c-3); the input bedrock wave error decreases from 0.113 to 0.081, and the waveform error at P1 decreases from 0.119 to 0.084; see Table 5). Finally, Fig. 3(b) shows the surface SI value distributions obtained when bC ref , C,(1) C,(2) best , and best are used as inputs. We can see that the 1st-order estimation gives SI distribution in high accuracy of 0.2/2.8 = 7%, and that the refined update improves the SI prediction to 3%. These comparisons demonstrate that ground motion distributions can be accurately estimated even for complex seismic waves with distinct properties (Kobe and Chuetsu-oki cases), and that the interactive refinement functions effectively for both cases.
18
4.4
Operational Timeline Analysis and Decision-Making Viability
Finally, we quantitatively evaluate the time-to-insight and computational resource savings from the moment a ground motion is observed at point P2 until decision-makers acquire the spatial time-history ground motion distribution for the Chuetsu-oki scenario. The actual post-disaster operational timeline progresses through the following steps (in practice, these steps are executed concurrently in parallel across numerous target locations): 1. Initial Surrogate Inference (1–2 s): Input observation data sC ref,P2 C,(1)
into the preconstructed N N P2 to infer the bedrock wave best . 2. 1st-Order Physical Simulation (12.1 min): Using the estimated bedrock wave, run a 3D nonlinear analysis on 128 nodes (128 GPUs) of Miyabi-G via the EBEGPU-Urgent kernel. A physically consistent firstorder spatial time-history distribution is obtained approximately 12 minutes post-event (completing immediate assessment). 3. On-the-fly Fine-Tuning (1.1 min): Based on the discrepancy between observed waves and the surface waves obtained from the 1st-order forward P2 analysis, construct N NC,(1) via dynamic incremental training on 1 node (1 GPU) of Miyabi-G. 4. Interactive Refined Simulation (12.1 min): Execute a second 3D C,(2) nonlinear analysis using the re-inferred bedrock wave best . Across the entire process above (1st-Order + Refinement), the pure computational time-to-insight is 25.3 minutes, and computational resource consumption is limited to 51.7 node-hours (128 nodes × 24.2 min + 1 node × 1.1 min) when refinement is triggered by automated thresholding. In conventional methods, achieving a comparable level of accuracy would require at least 102−3 forward analysis runs (approximately 20.2–202 hours on 128 nodes or 2,580–25,800 node-hours), even if utilizing the fast EBEGPU-Urgent kernel developed in this study, rendering real-time emergency operations physically impossible. The proposed workflow cuts required computation time and resource consumption by over 97.9% compared to conventional approaches while keeping physics-based simulations at its core. Consequently, it delivers 1st-order results just 12.1 minutes post-event while reliably satisfying the ”30-minute” urgent decisionmaking time frame, even when performing interactive refinement in response to decision-makers’ requirements.
5
Conclusion
In this study, we proposed a computational workflow for immediate post-disaster damage estimation that couples high-fidelity 3D nonlinear physics simulations with surrogate NNs, unifying normal-operation ”Resource Harvesting” with 19
emergency ”Urgent Response with Refinement.” Specifically, by developing CPUGPU-Harvest to leverage heterogeneous computing resources, we reduced energy-to-solution to 0.24 times that of conventional methods (a 76% reduction) and improved throughput 3.7-fold. This realized effective resource harvesting that flexibly utilizes short idle backfill slots during daily supercomputer operations, enabling the construction of high-precision pretraining datasets at a practical cost of 263 node-hours. Next, by coupling EBEGPU-Urgent—which prioritizes time-tosolution during emergencies—with surrogate NNs, we cut computation time and resource consumption by over 97.9% compared to conventional approaches. This paves the way to deliver physically consistent, scientifically sound insights— specifically high-fidelity spatial time-history ground motion distributions reflecting strong 3D nonlinear responses in large urban areas—within emergency decision-making timelines of 12.1 minutes post-event (1st-order estimation) and within 25.3 minutes (refined estimation) while incorporating interactive refinement. Future directions include further enhancing reliability by accounting for uncertainties in ground structures and more complex ground motion input conditions, as well as extending the framework to higher frequency ranges. Furthermore, beyond directly contributing to rapid, scientifically grounded decisionmaking for disaster mitigation during large-scale earthquakes, this study establishes an HPC utilization model bridging operational efficiency (”Harvesting”) and advanced decision-making (”Urgent Computing”) in future exa-scale environments, offering broad applicability to real-time inverse and identification analysis across diverse nonlinear time-evolution problems.
Appendix: High-fidelity physical model In the multi-spring model—one of the high-fidelity physical models used in 3D nonlinear ground amplification—incorporated in this study, physical properties are evaluated at four points per tetrahedral element using past history data to update the element stiffness matrix Ke composing Kn in Eqs. (1) P5 T T and (2) as Ke = j=1 wj Be,j De,j Be,j . Here, Be,j is a 6 × 30 matrix converting nodal displacements into strains, De,j is a 6 × 6 elastoplastic stiffness matrix at the integration point, and wj represents the weight of the integration point. In this P study, this elastoplastic stiffness matrix D is calcuN dτi ni nTi , where K is the bulk modulus and lated as D = KmmT + i=1 wis dγ i T m = {1, 1, 1, 0, 0, 0} . Additionally, γi , τi , and wis denote the strain, stress, and weight of the i-th 1D spring, respectively, and ni is a vector converting 1D spring strain into 3D strain. In other words, the multi-spring model is a physical model that evaluates 3D time-history nonlinear behavior by combining numerous experimentally derived 1D nonlinear springs, where strain induces nonlinearity, leading to stiffness degradation and increased damping. Here, because the modified Ramberg-Osgood model [27] and Masing’s rule [28] are employed as the physical model for each 1D spring, 40 bytes of data—comprising four double-precision variables and two flags—must be maintained per 1D spring. With the number of 1D springs set to N = 150, a total of 24 KB of history data 20
must be stored per tetrahedral element. Thus, while this enables high-fidelity physical simulations, maintaining such large history data results in a memory and computationally expensive simulation.
Appendix: Neural network for estimating input earthquake motion While various methods such as PINNs and FNO exist, considering stable performance guaranteed learning from a small dataset based on domain knowledge of earthquake engineering, we constructed a CNN-LSTM encoder-decoder network that estimates input bedrock waveforms (3-component x, y, z) from time-history surface observation waveforms (3-component x, y, z) recorded at observation points. As an architecture suited for handling both history dependence and dynamic locality/long-term temporal dependencies—widely used in waveform transformation tasks—we adopted a combination of CNN and LSTM [29]. Weight sharing in CNN enables efficient feature extraction in the time domain with fewer parameters, suppressing overfitting even when training data is limited [30]. First, target time-history waveforms were downsampled to the time domain corresponding to the target frequency range, yielding input and output time-history waveforms both of sequence length 1,813 with 3 components (1813 × 3 dimensions). The encoder transforms the input into latent features by stacking nc layers of 1D convolutional layers (kernel size k, stride 2) and ReLU activations, and extracts temporal dependencies using an LSTM (nLSTM layers, hidden dimension nhidden , same as the size of latent features). The decoder features a symmetric transposed convolutional structure; however, non-linear activation is omitted in the final layer, mapping the 3 components via grouped convolution (groups = 3), and the output length is matched to the input via linear interpolation. As described in the main text, N N P2 was pre-constructed using 100 pairs of random wave inputs and their corresponding responses at observation point P2 . Hyperparameters nc , nLSTM , nhidden , k, and the learning rate were optimized using Optuna [31] (nc ∈ {2, 3, 4}, nLSTM ∈ {1, 2, 3}, nhidden ∈ {128, 256, 512, 1024}, k ∈ {3, 5, 9, 17, 33, 65}, learning rate ∈ [5 × 10−5 , 5 × 10−4 ], up to 200 epochs per trial) using L1 loss and the Adam optimizer. Note that this was a singleobjective optimization targeting the validation loss, conducted over 100 trials using the TPE sampler. Using the 100 data pairs (80 training, 20 validation, 8 : 2 random split), early stopping was applied if the validation loss did not improve for 20 epochs, adopting the best-performing model. As a result of the optimization, nc = 2, nLSTM = 2, nhidden = 512, k = 3, and learning rate = 2.34 × 10−4 were obtained. The relative L1 loss on the validation data decreased to 0.0206 ± 0.0009 at an average of 159.6 epochs (hereafter, values represent the average across 50 random seeds, reported alongside 95% confidence intervals). Implemented in PyTorch [32], training was conducted on 1 node of NVIDIA GH200—envisioning pre-construction across multiple domains
21
P2 Figure 4: Combined training history of N N P2 and N NC,(1) (Chuetsu-oki scenario), separated by the dashed vertical line.
using supercomputer backfill resources—performing mini—batch training with a batch size of 10, requiring 97.6 minutes including optimization. To achieve event-specific adaptation, the pre-trained general-purpose model ∗,(1) ∗,(1) N N P2 is further trained on additional dataset (ssim,P2 and best ) to construct P2 N N∗,(1) . To prevent catastrophic forgetting—where overfitting to additional data destroys the broad representation capability acquired by the surrogate NN during pre-training [33] —Nrand = 20 random wave samples were mixed with Naug = 80 duplicated copies of the additional training sample, trained using mini-batches of size 10 and early stopping with patience=20 (evaluated on the loss of additional training samples). Fine-tuning for the Kobe and Chuetsu-oki waves required an average of 1.21 minutes (82.0 epochs) and 1.64 minutes (115.5 epochs), respectively. The achieved relative L1 errors were 0.0394 ± 0.0021 on the random wave validation data and 0.0328 ± 0.0012 on the additional training data for the Kobe wave, and 0.0307 ± 0.0019 on the random wave validation data and 0.0354 ± 0.0014 on the additional training data for the Chuetsu-oki wave. Figure 4 shows the concatenated training history of N N P2 and fine-tuning P2 history of N NC,(1) , with a gray dashed line marking the boundary between the two (Chuetsu-oki scenario). The blue line represents the error on the random wave validation data across both sides of the dashed line (left: N N P2 , right: P2 N NC,(1) ), while the orange line represents the error on the additional data for P2 N NC,(1) on the right side of the dashed line. While the fitting performance (orange line) improves through fine-tuning, the random wave validation error (blue line) degrades only slightly, confirming that generalization performance is not catastrophically impaired. It should be noted that the training/validation C,(1) C,(1) data consisted solely of random waves and (ssim,P2 , best ), and the test data (true values) was not used during training or hyperparameter tuning.
22
Acknowledgment The authors thank the Mainline Committee of the Association for the Development of Earthquake Prediction (ADEP) for providing the ground model used in this work. This work was supported by JSPS KAKENHI (Grant Numbers 26H02178, 25K21686).
References [1] M. D. Trifunac, ”Site conditions and earthquake ground motion – A review,” Soil Dynamics and Earthquake Engineering, vol. 90, pp. 88–100, 2016. https://doi.org/10.1016/j.soildyn.2016.08.003 [2] J. Liang and S. Sun, ”Site effects on seismic behavior of pipelines: A review,” ASME Journal of Pressure Vessel Technology, vol. 122, no. 4, pp. 469–475, 2000. https://doi.org/10.1115/1.1285974 [3] U.S. Geological Survey (USGS), ”ShakeMap,” https://earthquake. usgs.gov/data/shakemap/ (accessed Aug. 10, 2026). [4] National Research Institute for Earth Science and Disaster Resilience (NIED), ”J-RISQ: Real-time Earthquake Information System for Disaster Response,” https://www.j-risq.bosai.go.jp/report/en/ (accessed Aug. 10, 2026). [5] Caltech and UCLA, ”Community Seismic Network,” http://csn. caltech.edu/ (accessed Aug. 10, 2026). [6] B. Enders, D. Bard, C. Snavely, L. Gerhardt, J. Lee, B. Totzke, K. Antypas, S. Byna, R. Cheema, S. Cholia, M. Day, A. Gaur, A. Greiner, T. Groves, M. Kiran, Q. Koziol, K. Rowland, C. Samuel, A. Selvarajan, A. Sim, D. Skinner, R. Thomas, and G. Torok, ”Cross-facility science with the Superfacility Project at LBNL,” in Proc. IEEE/ACM 2nd Annual Workshop on Extreme-scale Experiment-in-the-Loop Computing (XLOOP), pp. 1–7, 2020. https://doi.org/10.1109/XLOOP51963.2020.00006 [7] T. Goubier, N. Rakowsky, and S. Harig, ”Fast tsunami simulations for a real-time emergency response flow,” in Proc. IEEE/ACM HPC for Urgent Decision Making (UrgentHPC), pp. 21–26, 2020. https://doi.org/10. 1109/UrgentHPC51945.2020.00008 [8] S. Henneking, S. Venkat, V. Dobrev, J. Camier, T. Kolev, M. Fernando, A.-A. Gabriel, and O. Ghattas, ”Real-time Bayesian inference at extreme scale: A digital twin for tsunami early warning applied to the Cascadia Subduction Zone,” in Proc. Int. Conf. High Performance Computing, Networking, Storage and Analysis (SC), pp. 60–71, 2025. https: //doi.org/10.1145/3712285.3771787
23
[9] T. Ichimura, K. Fujita, S. Tanaka, M. Hori, M. Lalith, Y. Shizawa, and H. Kobayashi, ”Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured FE non-linear seismic wave simulation,” in Proc. Int. Conf. High Performance Computing, Networking, Storage and Analysis (SC), pp. 15–26, 2014. https://doi.org/10.1109/ SC.2014.7 [10] T. Ichimura, K. Fujita, P. E. B. Quinay, M. Lalith, M. Hori, S. Tanaka, Y. Shizawa, H. Kobayashi, and K. Minami, ”Implicit nonlinear wave simulation with 1.08T DOF and 0.270T unstructured finite elements to enhance comprehensive earthquake simulation,” in Proc. Int. Conf. High Performance Computing, Networking, Storage and Analysis (SC), pp. 1–12, 2015. https://doi.org/10.1145/2807591.2807674 [11] T. Ichimura, K. Fujita, M. Hori, T. Sakanoue, and R. Hamanaka, ”Threedimensional nonlinear seismic ground response analysis of local site effects for estimating seismic behavior of buried pipelines,” ASME Journal of Pressure Vessel Technology, vol. 136, no. 4, Art. no. 041702, 2014. https://doi.org/10.1115/1.4026208 [12] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics, 6th ed. Amsterdam, The Netherlands: Elsevier, 2005. [13] N. M. Newmark, ”A method of computation for structural dynamics,” Journal of the Engineering Mechanics Division, ASCE, vol. 85, no. 3, pp. 67–94, 1959. https://doi.org/10.1061/JMCEA3.0000098 [14] S. Iai, ”Three dimensional formulation and objectivity of a strain space multiple mechanism model for sand,” Soils and Foundations, vol. 33, no. 1, pp. 192–199, 1993. https://doi.org/10.3208/sandf1972.33.192 [15] T. Ichimura, K. Fujita, H. Ito, M. Hori, and L. Maddegedara, ”Accelerating nonlinear time-history analysis with complex constitutive laws via heterogeneous memory management: From 3D seismic simulation to neural network training,” in Lecture Notes in Computer Science, vol. 16784, Springer, pp. 257–271, 2026. https://doi.org/10.1007/978-3-032-29924-6_19 [16] PCI-SIG, ”PCIe 5.0 Specification,” https://pcisig.com/ (accessed Aug. 10, 2026). [17] NVIDIA, ”NVIDIA GH200 Grace Hopper Superchip Architecture,” https://nvdam.widen.net/s/c9lts6msjj/ nvidia-grace-hopper-superchip-architecture-whitepaper (accessed Aug. 10, 2026). [18] J. M. Winget and T. J. R. Hughes, ”Solution algorithms for nonlinear transient heat conduction analysis employing element-by-element iterative strategies,” Computer Methods in Applied Mechanics and Engi24
neering, vol. 52, no. 1–3, pp. 711–815, 1985. https://doi.org/10.1016/ 0045-7825(85)90011-3 [19] C. Augonnet, S. Thibault, R. Namyst, and P.-A. Wacrenier, ”StarPU: A unified platform for task scheduling on heterogeneous multicore architectures,” Concurrency and Computation: Practice and Experience, vol. 23, no. 2, pp. 187–198, 2011. https://doi.org/10.1002/cpe.1631 [20] O. Pearce, ”Experiences using CPUs and GPUs for cooperative computation in a multi-physics simulation,” in Proc. 47th Int. Conf. on Parallel Processing Workshops (ICPP Workshops), Art. no. 25, pp. 1–10, 2018. https://doi.org/10.1145/3229710.3229711 [21] S. Mittal and J. S. Vetter, ”A survey of CPU-GPU heterogeneous computing techniques,” ACM Computing Surveys, vol. 47, no. 4, Art. no. 69, 2015. https://doi.org/10.1145/2788396 [22] K. Raju and N. N. Chiplunkar, ”A survey on techniques for cooperative CPU-GPU computing,” Sustainable Computing: Informatics and Systems, vol. 19, pp. 72–85, 2018. https://doi.org/10.1016/j.suscom.2018.07. 010 [23] Information Technology Center, The University of Tokyo, ”Supercomputer System ‘Miyabi’,” https://www.cc.u-tokyo.ac.jp/en/supercomputer/ miyabi/service/ (accessed Aug. 10, 2026). [24] T. Ichimura, K. Fujita, M. Hori, L. Maddegedara, J. Wells, A. Gray, I. Karlin, and J. Linford, ”Heterogeneous computing in a strongly-connected CPU-GPU environment: Fast multiple time-evolution equation-based modeling accelerated using data-driven approach,” in Proc. IEEE/ACM Int. Conf. for High Performance Computing, Networking, Storage and Analysis Workshops (SC Workshops), pp. 1967–1978, 2024. https://doi.org/10. 1109/SCW63240.2024.00246 [25] Japan Meteorological Agency, ”Kyoshin Seismic Data,” https://www. data.jma.go.jp/eqev/data/kyoshin/jishin/ (accessed Aug. 10, 2026). [26] G. W. Housner, ”Spectrum intensities of strong-motion earthquakes,” in Proc. Symposium on Earthquakes and Blast Effects on Structures, pp. 20– 36, 1952. [27] I. M. Idriss, R. Dobry, and R. D. Singh, ”Nonlinear behavior of soft clays during cyclic loading,” Journal of the Geotechnical Engineering Division, ASCE, vol. 104, no. 12, pp. 1427–1447, 1978. https://doi.org/10.1061/ AJGEB6.0000727 [28] G. Masing, ”Eigenspannungen und Verfestigung beim Messing,” in Proc. 2nd Int. Congr. of Applied Mechanics, pp. 332–335, 1926.
25
[29] K. Tan, D. Wang, ”A Convolutional Recurrent Neural Network for RealTime Speech Enhancement”, Interspeech 2018. https://doi:10.21437/ Interspeech.2018-1405 [30] Y. Lecun, L. Bottou, Y. Bengio, P. Haffner, ”Gradient-based learning applied to document recognition,” Proc. of the IEEE 86.11 pp. 2278–2324, 1998. https://doi.org/10.1109/5.726791 [31] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, ”Optuna: A nextgeneration hyperparameter optimization framework,” in Proc. 25th ACM SIGKDD Int. Conf. on Knowledge Discovery & Data Mining (KDD), pp. 2623–2631, 2019. https://doi.org/10.1145/3292500.3330701 [32] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, ”PyTorch: An imperative style, high-performance deep learning library,” in Advances in Neural Information Processing Systems (NeurIPS), vol. 32, pp. 8026–8037, 2019. [33] A. Robins, ”Catastrophic forgetting, rehearsal and pseudorehearsal,” Connection Science 7.2 pp 123–146, 1995. https://doi.org/10.1080/ 09540099550039318
26