From 57472f8e428c1f8b3e53412b90ca5999787d69f8 Mon Sep 17 00:00:00 2001 From: Andrew Hicks Date: Sat, 29 Apr 2023 18:54:14 -0500 Subject: [PATCH] Cleaned up qtensor3d --- q3d/compute.py | 85 +++++++++++++++++++++++------------------------ q3d/config.py | 4 ++- scripts/qtensor3d | 85 ++++++++++++++++++++++++++--------------------- 3 files changed, 92 insertions(+), 82 deletions(-) diff --git a/q3d/compute.py b/q3d/compute.py index a793412..04e431e 100644 --- a/q3d/compute.py +++ b/q3d/compute.py @@ -1,47 +1,6 @@ from sympyplus import * -from q3d.config import constants as c -from q3d.config import settings - -# Set up Qvector objects - -nu = AbstractVector('nu') - -q = QVector('q') -Dq = q.grad -Dqq = Param([Dq,q]) -Q = q.tens - -p = QVector('p') -Dp = p.grad -Dpp = Param([Dp,p]) -P = p.tens - -r = QVector('r') -Dr = r.grad -Drr = Param([Dr,r]) -R = r.tens - -qp = QVector('q_prev') -Dqp = qp.grad -QP = qp.tens - -qpp = QVector('q_prev_prev') -Dqpp = qpp.grad -QPP = qpp.tens - -qnp = QVector('q_newt_prev') -Dqnp = qnp.grad -Dqnpqnp = Param([Dqnp,qnp]) -QNP = qnp.tens - -f = QVector('f') -g = QVector('g') - -Q0 = c.S0*(outerp(nu,nu) - (1.0/3.0)*eye(3)) -Pi = eye(3) - outerp(nu,nu) - -################################################################# +# functions that give the individual terms def termL1(grad1,grad2): return innerp(grad1,grad2) @@ -71,9 +30,49 @@ def term_twist_var(q,p): return 2*c.q0*mixedp(Q,P) + 2*c.q0*mixedp(P,Q) -######################################################### +# compute function, which puts the PDE system together def compute(): + from q3d.config import constants as c + from q3d.config import settings + + # set up Qvector objects + nu = AbstractVector('nu') + + q = QVector('q') + Dq = q.grad + Dqq = Param([Dq,q]) + Q = q.tens + + p = QVector('p') + Dp = p.grad + Dpp = Param([Dp,p]) + P = p.tens + + r = QVector('r') + Dr = r.grad + Drr = Param([Dr,r]) + R = r.tens + + qp = QVector('q_prev') + Dqp = qp.grad + QP = qp.tens + + qpp = QVector('q_prev_prev') + Dqpp = qpp.grad + QPP = qpp.tens + + qnp = QVector('q_newt_prev') + Dqnp = qnp.grad + Dqnpqnp = Param([Dqnp,qnp]) + QNP = qnp.tens + + f = QVector('f') + g = QVector('g') + + Q0 = c.S0*(outerp(nu,nu) - (1.0/3.0)*eye(3)) + Pi = eye(3) - outerp(nu,nu) + # Define 'energies' used to calculate the energy energies = EnergyForm(Dqq,Dpp,Drr) energies.add_domain(GeneralForm(c.L1/2*termL1(Dq,Dq)+c.L2/2*termL2(Dq,Dq)+c.L3/2*termL3(Dq,Dq),Dqq,name='Elastic Energy')) diff --git a/q3d/config.py b/q3d/config.py index 02e714c..c788266 100644 --- a/q3d/config.py +++ b/q3d/config.py @@ -67,7 +67,7 @@ def nondimensionalize(c:FromDict,R:float): except AttributeError: c.beta = 1 -def initialize(settings_path, constants_path=None, supersessions={}): +def initialize(settings_path, constants_path=None, *, supersessions={}): from q3d.loaddump import load_yml # make settings global, thus importable @@ -104,6 +104,8 @@ def initialize(settings_path, constants_path=None, supersessions={}): # process constants process_constants(constants) + return settings, constants + def main(): initialize('settings/settings.yml','constants/5cb.yml') print(constants) diff --git a/scripts/qtensor3d b/scripts/qtensor3d index 5818ffd..678da48 100644 --- a/scripts/qtensor3d +++ b/scripts/qtensor3d @@ -3,12 +3,20 @@ import getopt import os import sys -from q3d.loaddump import load_json -import q3d.uflcache as uflcache +from time import sleep + from firedrake import COMM_WORLD as comm -from q3d.firedrakeplus import ManuQ -from q3d.firedrakeplus import set_eqn_globals, solve_PDE, compute_energy, errorL2, errorH1, choose_mesh + +import q3d.compute as compute +import q3d.config as config import q3d.printoff as pr +import q3d.saves as saves +import q3d.uflcache as uflcache +from q3d.firedrakeplus import (ManuQ, choose_mesh, compute_energy, errorH1, + errorL2, set_eqn_globals, solve_PDE) +from q3d.loaddump import load_json +from q3d.misc import Timer, check, get_range + def print0(*args,**kwargs): if comm.rank == 0: @@ -52,19 +60,32 @@ def check_if_checkpoint_exists(path): print0("cannot resume since no checkpoint found. Try overwriting instead.") sys.exit() -def run(path, *, mode='r', supersessions={}): - # these three modules must be imported in order and before other modules, or else they won't work properly - import q3d.config as config - config.initialize(f'{path}/settings.yml',f'{path}/constants.yml',supersessions) - from q3d.config import settings - from time import sleep +def run(path, *, overwrite=False, supersessions={}): + # first, check if path is a valid save + check_if_valid_save(path) - import q3d.saves as saves + # then, check if checkpoint exists if in resume mode + if not overwrite: + check_if_checkpoint_exists(path) + + # if in overwrite mode, make sure the user really wants to overwrite + if overwrite and not answers_yes_to(f"will overwrite save at '{path}'. Are you sure you want to continue? (y/n) "): + pr.text('exiting') + return + + # confirm th action of overwriting or resuming + BLANK = 'overwrit' if overwrite else 'resum' + pr.info(f"{BLANK}ing save at '{path}'") + + # must initialize q3d.config before q3d.saves, otherwise will break + settings, constants = config.initialize(f'{path}/settings.yml', f'{path}/constants.yml', supersessions=supersessions) + mode = 'o' if overwrite else 'r' saves.initialize(mode,path) - # import other modules - import q3d.printoff as pr - from q3d.misc import Timer, get_range, check + # if we have time steps that, when added together, would be 16 digits or more, we must stop this + if (dt := settings.time.step) < 1.0e-5 or dt >= 1.0e+9: + pr.fail(f'step size {dt} not allowed') + return # print info pr.constants_info() @@ -79,14 +100,11 @@ def run(path, *, mode='r', supersessions={}): # perform sympyplus preliminary computations timer = Timer() timer.start() - - import q3d.compute as compute comp = compute.compute() - timer.stop() # rebuild UFL cache if in overwrite mode - if mode == 'o': + if overwrite: pr.text("Rebuilding UFL cache...",end=' ') uflcache.build_uflcache(path) pr.text("build successful.") @@ -176,27 +194,22 @@ def auto_run(path, *, overwrite=False, supersessions): supersessions['save-every'] = supersessions.get('save-every', 50) if 'no-gd' in supersessions: print0("cannot choose --no-gd in auto mode") - sys.exit() + return for i in range(100): + # assume not completed, then run iterations until completed, adjusting the step size as appropriate completed = False while not completed: - if i == 0 and overwrite: - completed = overwrite_save(path, supersessions=supersessions) - else: - completed = resume_save(path, supersessions=supersessions) - - if completed: - supersessions['dt'] *= 10 if supersessions['dt'] <= 1_000_000_000 else 1 - continue - - supersessions['dt'] *= 0.1 + completed = run(path, overwrite=overwrite, supersessions=supersessions) + supersessions['dt'] *= 10 if completed else 0.1 + # turn off overwrite mode after iteration i == 0 is completed + overwrite = False def main(): sys_argv = sys.argv[1:] if len(sys_argv) == 0: print0('no argument supplied') - sys.exit() + return try: opts, listargs = getopt.getopt(sys.argv[1:],'oa',['no-gd','dt=','num-steps=','save-every=','help']) @@ -204,7 +217,7 @@ def main(): # print help information and exit: print0(err) print0('use --help for usage') - sys.exit() + return # resume mode is default overwrite = False @@ -225,20 +238,16 @@ def main(): # get help elif o in ('--help'): usage() - sys.exit() + return else: - sys.exit() + return # overwrite or resume if auto: auto_run(listargs[-1], overwrite=overwrite, supersessions=supersessions) return - if overwrite: - overwrite_save(listargs[-1], supersessions=supersessions) - return - - resume_save(listargs[-1], supersessions=supersessions) + run(listargs[-1], overwrite=overwrite, supersessions=supersessions) if __name__ == '__main__': main()