From 3523d1da0b672e72262b0857797a4c31f1179ba3 Mon Sep 17 00:00:00 2001 From: Vatsal Sanjay Date: Sun, 3 Nov 2024 13:17:00 +0100 Subject: [PATCH] =?UTF-8?q?=F0=9F=92=8E=20feat(viscoelastic-scalar-3D):=20?= =?UTF-8?q?Refactor=20tensor=20operations=20and=20improve=20maintainabilit?= =?UTF-8?q?y?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Refactored tensor operations to use an intermediate tensor structure `A` for improved clarity and maintainability - Unified the `A`-tensor-based approach throughout the code for consistent tensor manipulation - Added explicit symmetry enforcement in tensor operations - Improved readability by separating tensor transformation steps - Enhanced extensibility for future tensor-only implementations --- .../log-conform-viscoelastic-scalar-3D.h | 67 +++++++++++++------ 1 file changed, 46 insertions(+), 21 deletions(-) diff --git a/src-local/log-conform-viscoelastic-scalar-3D.h b/src-local/log-conform-viscoelastic-scalar-3D.h index dd5c25a..a353028 100644 --- a/src-local/log-conform-viscoelastic-scalar-3D.h +++ b/src-local/log-conform-viscoelastic-scalar-3D.h @@ -1,12 +1,12 @@ /** Title: log-conform-viscoelastic-3D.h -# Version: 2.1 +# Version: 2.2 # Main feature 1: A exists in across the domain and relaxes according to \lambda. The stress only acts according to G. # Main feature 2: This is the 3D implementation of [log-conform-viscoelastic-scalar-2D.h](log-conform-viscoelastic-scalar-2D.h). # Author: Vatsal Sanjay # vatsalsanjay@gmail.com # Physics of Fluids -# Updated: Oct 29, 2024 +# Updated: Nov 3, 2024 # change log: Oct 19, 2024 (v1.0) - 3D implementation @@ -30,6 +30,13 @@ # change log: Oct 29, 2024 (v2.1) - Added some initialization functions for pseudo_v and pseudo_t and their 3D counterparts. + +# change log: Nov 3, 2024 (v2.2) +- Refactored tensor operations to use intermediate tensor structure A for improved clarity and maintainability +- Unified A-tensor-based approach throughout the code for consistent tensor manipulation +- Added explicit symmetry enforcement in tensor operations +- Improved readability by separating tensor transformation steps +- Enhanced extensibility for future tensor-only implementations */ /** The code is same as http://basilisk.fr/src/log-conform.h but @@ -848,30 +855,48 @@ event tracer_advection(i++) Lambda.z = exp(Lambda.z); // Reconstruct A using A = R * diag(Lambda) * R^T - A11[] = Lambda.x * sq(R.x.x) + Lambda.y * sq(R.x.y) + Lambda.z * sq(R.x.z); - A12[] = Lambda.x * R.x.x * R.y.x + Lambda.y * R.x.y * R.y.y + Lambda.z * R.x.z * R.y.z; - A13[] = Lambda.x * R.x.x * R.z.x + Lambda.y * R.x.y * R.z.y + Lambda.z * R.x.z * R.z.z; - A22[] = Lambda.x * sq(R.y.x) + Lambda.y * sq(R.y.y) + Lambda.z * sq(R.y.z); - A23[] = Lambda.x * R.y.x * R.z.x + Lambda.y * R.y.y * R.z.y + Lambda.z * R.y.z * R.z.z; - A33[] = Lambda.x * sq(R.z.x) + Lambda.y * sq(R.z.y) + Lambda.z * sq(R.z.z); + A.x.x = Lambda.x * sq(R.x.x) + Lambda.y * sq(R.x.y) + Lambda.z * sq(R.x.z); + A.x.y = Lambda.x * R.x.x * R.y.x + Lambda.y * R.x.y * R.y.y + Lambda.z * R.x.z * R.y.z; + A.y.x = A.x.y; + A.x.z = Lambda.x * R.x.x * R.z.x + Lambda.y * R.x.y * R.z.y + Lambda.z * R.x.z * R.z.z; + A.z.x = A.x.z; + A.y.y = Lambda.x * sq(R.y.x) + Lambda.y * sq(R.y.y) + Lambda.z * sq(R.y.z); + A.y.z = Lambda.x * R.y.x * R.z.x + Lambda.y * R.y.y * R.z.y + Lambda.z * R.y.z * R.z.z; + A.z.y = A.y.z; + A.z.z = Lambda.x * sq(R.z.x) + Lambda.y * sq(R.z.y) + Lambda.z * sq(R.z.z); // Apply relaxation using the relaxation time lambda double intFactor = lambda[] != 0. ? exp(-dt/lambda[]) : 0.; - - A11[] = 1. + (A11[] - 1.)*intFactor; - A22[] = 1. + (A22[] - 1.)*intFactor; - A33[] = 1. + (A33[] - 1.)*intFactor; - A12[] *= intFactor; - A13[] *= intFactor; - A23[] *= intFactor; + + A.x.y *= intFactor; + A.y.x = A.x.y; + A.x.z *= intFactor; + A.z.x = A.x.z; + A.y.z *= intFactor; + A.z.y = A.y.z; + foreach_dimension() + A.x.x = 1. + (A.x.x - 1.)*intFactor; + + /* + Get Aij from A. These commands might look repetitive. But, I do this so that in the future, generalization to tensor only form is easier. + */ + + // diagonal terms: + A11[] = A.x.x; + A22[] = A.y.y; + A33[] = A.z.z; + // off-diagonal terms: + A12[] = A.x.y; + A13[] = A.x.z; + A23[] = A.y.z; // Compute the stress tensor T using the polymer modulus Gp - T11[] = Gp[]*(A11[] - 1.); - T22[] = Gp[]*(A22[] - 1.); - T33[] = Gp[]*(A33[] - 1.); - T12[] = Gp[]*A12[]; - T13[] = Gp[]*A13[]; - T23[] = Gp[]*A23[]; + T11[] = Gp[]*(A.x.x - 1.); + T22[] = Gp[]*(A.y.y - 1.); + T33[] = Gp[]*(A.z.z - 1.); + T12[] = Gp[]*A.x.y; + T13[] = Gp[]*A.x.z; + T23[] = Gp[]*A.y.z; } } #endif